Move the csc directory into spm, because it handles all csc, ijv, csr formats

CMakeLists.txt

+# include(RulesPrecisions)
+# # reset variables
+# set(generated_files "")
+# set(generated_headers "")
+
+# set(HEADERS
+#   z_spm.h
+# )
+
+# ### generate the dsparse_cores headers for all possible precisions
+# precisions_rules_py(generated_headers
+#                     "${HEADERS}"
+#                     PRECISIONS "p;s;d;c;z")
+
+# add_custom_target(spm_headers
+#   DEPENDS ${generated_headers} )
+

csc.h

+/**
+ *
+ * @file csc.h
+ *
+ *  PaStiX sparse matrix routines to handle different format of sparse matrices.
+ *  $COPYRIGHTS$
+ *
+ * @version 5.1.0
+ * @author Xavier Lacoste
+ * @author Pierre Ramet
+ * @author Mathieu Faverge
+ * @date 2013-06-24
+ *
+ **/
+#ifndef _CSC_H_
+#define _CSC_H_
+
+/**
+ * @ingroup pastix_spm
+ *
+ * @struct pastix_spm_s - Sparse matrix data structure
+ */
+struct pastix_spm_s {
+    int           mtxtype;   /*< Matrix structure: PastixGeneral, PastixSymmetric or PastixHermitian.          */
+    int           flttype;   /*< avals datatype: PastixPattern, PastixFloat, PastixDouble, PastixComplex32 or PastixComplex64 */
+    int           fmttype;   /*< Matrix storage format: PastixCSC, PastixCSR, PastixIJV                        */
+    pastix_int_t  gN;        /*< Global number of vertices in the compressed graph   */
+    pastix_int_t  n;         /*< Local number of vertices in the compressed graph    */
+    pastix_int_t  gnnz;      /*< Global number of non zeroes in the compressed graph */
+    pastix_int_t  nnz;       /*< Local number of non zeroes in the compressed graph  */
+    pastix_int_t  dof;       /*< Number of degrees of freedom per unknown            */
+    pastix_int_t *colptr;    /*< List of indirections to rows for each vertex        */
+    pastix_int_t *rowptr;    /*< List of edges for each vertex                       */
+    pastix_int_t *loc2glob;  /*< Corresponding numbering from local to global        */
+    void         *values;    /*< Values stored in the matrix                         */
+};
+typedef struct pastix_spm_s pastix_csc_t;
+
+int
+csc_load( pastix_int_t  *n,
+          pastix_int_t **colptr,
+          pastix_int_t **rows,
+          int           *valtype,
+          void         **values,
+          int           *dof,
+          FILE          *infile );
+
+int
+csc_save( pastix_int_t  n,
+          pastix_int_t *colptr,
+          pastix_int_t *rows,
+          int           ft,
+          void         *values,
+          int           dof,
+          FILE         *outfile );
+
+int cscLoad( pastix_csc_t *csc, FILE *infile );
+int cscSave( pastix_csc_t *csc, FILE *outfile );
+
+
+int spmGenRHS(int type, int nrhs, const pastix_csc_t *spm, void *x, int ldx, void *b, int ldb );
+int spmCheckAxb( int nrhs, const pastix_csc_t *spm, void *x0, int ldx0, void *b, int ldb, const void *x, int ldx );
+
+void          spmInit( pastix_csc_t *spm );
+void          spmExit( pastix_csc_t *spm );
+pastix_csc_t *spmCopy( const pastix_csc_t *spm );
+void          spmBase( pastix_csc_t *spm, int baseval );
+int           spmConvert( int ofmttype, pastix_csc_t *ospm );
+pastix_int_t  spmFindBase( const pastix_csc_t *spm );
+double        spmNorm( int ntype, const pastix_csc_t *csc );
+int           spmMatVec(int trans, const void *alpha, const pastix_csc_t *csc, const void *x, const void *beta, void *y );
+
+int           spmSort( pastix_csc_t *csc );
+pastix_int_t  spmMergeDuplicate( pastix_csc_t *csc );
+pastix_int_t  spmSymmetrize( pastix_csc_t *csc );
+
+pastix_csc_t *spmCheckAndCorrect( pastix_csc_t *csc );
+
+#endif /* _CSC_H_ */

csc_io.c

+/**
+ *
+ * @file csc_io.c
+ *
+ *  PaStiX csc routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Xavier Lacoste
+ * @author Pierre Ramet
+ * @author Mathieu Faverge
+ * @date 2013-06-24
+ *
+ **/
+#include "common.h"
+#include "csc.h"
+
+static inline int
+readArrayOfInteger( FILE         *stream,
+                    pastix_int_t  n,
+                    pastix_int_t *array )
+{
+    long tmp1, tmp2, tmp3, tmp4;
+    pastix_int_t i;
+
+    /* Read 4 by 4 */
+    for (i=0; i<(n-3); i+=4)
+    {
+        if (4 != fscanf(stream, "%ld %ld %ld %ld", &tmp1, &tmp2, &tmp3, &tmp4)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+
+        array[i  ] = (pastix_int_t)tmp1;
+        array[i+1] = (pastix_int_t)tmp2;
+        array[i+2] = (pastix_int_t)tmp3;
+        array[i+3] = (pastix_int_t)tmp4;
+    }
+
+    assert( n-i < 4 );
+    switch ( n - i )
+    {
+    case 3:
+        if (3 != fscanf(stream, "%ld %ld %ld", &tmp1, &tmp2, &tmp3)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+
+        array[i  ] = (pastix_int_t)tmp1;
+        array[i+1] = (pastix_int_t)tmp2;
+        array[i+2] = (pastix_int_t)tmp3;
+        break;
+    case 2:
+        if (2 != fscanf(stream, "%ld %ld", &tmp1, &tmp2)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+
+        array[i  ] = (pastix_int_t)tmp1;
+        array[i+1] = (pastix_int_t)tmp2;
+        break;
+    case 1:
+        if (1 != fscanf(stream, "%ld", &tmp1)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+
+        array[i  ] = (pastix_int_t)tmp1;
+        break;
+    }
+
+    return EXIT_SUCCESS;
+}
+
+static inline int
+readArrayOfComplex64( FILE               *stream,
+                      pastix_int_t        n,
+                      pastix_complex64_t *array )
+{
+    double tmp1, tmp2, tmp3, tmp4;
+    double tmp5, tmp6, tmp7, tmp8;
+    pastix_int_t i;
+
+    /* Read 4 by 4 */
+    for (i=0; i<(n-3); i+=4)
+    {
+        if (8 != fscanf(stream, "%lg %lg %lg %lg %lg %lg %lg %lg",
+                        &tmp1, &tmp2, &tmp3, &tmp4,
+                        &tmp5, &tmp6, &tmp7, &tmp8)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (pastix_complex64_t)(tmp1 + I * tmp2);
+        array[i+1] = (pastix_complex64_t)(tmp3 + I * tmp4);
+        array[i+2] = (pastix_complex64_t)(tmp5 + I * tmp6);
+        array[i+3] = (pastix_complex64_t)(tmp7 + I * tmp8);
+    }
+
+    assert( n-i < 4 );
+    switch ( n - i )
+    {
+    case 3:
+        if (6 != fscanf(stream, "%lg %lg %lg %lg %lg %lg",
+                        &tmp1, &tmp2, &tmp3, &tmp4,
+                        &tmp5, &tmp6)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (pastix_complex64_t)(tmp1 + I * tmp2);
+        array[i+1] = (pastix_complex64_t)(tmp3 + I * tmp4);
+        array[i+2] = (pastix_complex64_t)(tmp5 + I * tmp6);
+        break;
+
+    case 2:
+        if (4 != fscanf(stream, "%lg %lg %lg %lg",
+                        &tmp1, &tmp2, &tmp3, &tmp4)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (pastix_complex64_t)(tmp1 + I * tmp2);
+        array[i+1] = (pastix_complex64_t)(tmp3 + I * tmp4);
+        break;
+
+    case 1:
+        if (2 != fscanf(stream, "%lg %lg",
+                        &tmp1, &tmp2)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (pastix_complex64_t)(tmp1 + I * tmp2);
+        break;
+    }
+
+    return EXIT_SUCCESS;
+}
+
+static inline int
+readArrayOfComplex32( FILE               *stream,
+                      pastix_int_t        n,
+                      pastix_complex32_t *array )
+{
+    float tmp1, tmp2, tmp3, tmp4;
+    float tmp5, tmp6, tmp7, tmp8;
+    pastix_int_t i;
+
+    /* Read 4 by 4 */
+    for (i=0; i<(n-3); i+=4)
+    {
+        if (8 != fscanf(stream, "%g %g %g %g %g %g %g %g",
+                        &tmp1, &tmp2, &tmp3, &tmp4,
+                        &tmp5, &tmp6, &tmp7, &tmp8)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (pastix_complex32_t)(tmp1 + I * tmp2);
+        array[i+1] = (pastix_complex32_t)(tmp3 + I * tmp4);
+        array[i+2] = (pastix_complex32_t)(tmp5 + I * tmp6);
+        array[i+3] = (pastix_complex32_t)(tmp7 + I * tmp8);
+    }
+
+    assert( n-i < 4 );
+    switch ( n - i )
+    {
+    case 3:
+        if (6 != fscanf(stream, "%g %g %g %g %g %g",
+                        &tmp1, &tmp2, &tmp3, &tmp4,
+                        &tmp5, &tmp6)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (pastix_complex32_t)(tmp1 + I * tmp2);
+        array[i+1] = (pastix_complex32_t)(tmp3 + I * tmp4);
+        array[i+2] = (pastix_complex32_t)(tmp5 + I * tmp6);
+        break;
+
+    case 2:
+        if (4 != fscanf(stream, "%g %g %g %g",
+                        &tmp1, &tmp2, &tmp3, &tmp4)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (pastix_complex32_t)(tmp1 + I * tmp2);
+        array[i+1] = (pastix_complex32_t)(tmp3 + I * tmp4);
+        break;
+
+    case 1:
+        if (2 != fscanf(stream, "%g %g",
+                        &tmp1, &tmp2)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (pastix_complex32_t)(tmp1 + I * tmp2);
+        break;
+    }
+
+    return EXIT_SUCCESS;
+}
+
+static inline int
+readArrayOfDouble( FILE         *stream,
+                   pastix_int_t  n,
+                   double       *array )
+{
+    double tmp1, tmp2, tmp3, tmp4;
+    pastix_int_t i;
+
+    /* Read 4 by 4 */
+    for (i=0; i<(n-3); i+=4)
+    {
+        if (4 != fscanf(stream, "%lg %lg %lg %lg",
+                        &tmp1, &tmp2, &tmp3, &tmp4)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (double)(tmp1);
+        array[i+1] = (double)(tmp2);
+        array[i+2] = (double)(tmp3);
+        array[i+3] = (double)(tmp4);
+    }
+
+    assert( n-i < 4 );
+    switch ( n - i )
+    {
+    case 3:
+        if (1 != fscanf(stream, "%lg %lg %lg",
+                        &tmp1, &tmp2, &tmp3)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (double)(tmp1);
+        array[i+1] = (double)(tmp2);
+        array[i+2] = (double)(tmp3);
+        break;
+
+    case 2:
+        if (2 != fscanf(stream, "%lg %lg",
+                        &tmp1, &tmp2)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (double)(tmp1);
+        array[i+1] = (double)(tmp2);
+        break;
+
+    case 1:
+        if (1 != fscanf(stream, "%lg",
+                        &tmp1)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (double)(tmp1);
+        break;
+    }
+
+    return EXIT_SUCCESS;
+}
+
+static inline int
+readArrayOfFloat( FILE         *stream,
+                  pastix_int_t  n,
+                  float        *array )
+{
+    float tmp1, tmp2, tmp3, tmp4;
+    pastix_int_t i;
+
+    /* Read 4 by 4 */
+    for (i=0; i<(n-3); i+=4)
+    {
+        if (4 != fscanf(stream, "%g %g %g %g",
+                        &tmp1, &tmp2, &tmp3, &tmp4)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (float)(tmp1);
+        array[i+1] = (float)(tmp2);
+        array[i+2] = (float)(tmp3);
+        array[i+3] = (float)(tmp4);
+    }
+
+    assert( n-i < 4 );
+    switch ( n - i )
+    {
+    case 3:
+        if (3 != fscanf(stream, "%g %g %g",
+                        &tmp1, &tmp2, &tmp3)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (float)(tmp1);
+        array[i+1] = (float)(tmp2);
+        array[i+2] = (float)(tmp3);
+        break;
+
+    case 2:
+        if (2 != fscanf(stream, "%g %g",
+                        &tmp1, &tmp2)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (float)(tmp1);
+        array[i+1] = (float)(tmp2);
+        break;
+
+    case 1:
+        if (1 != fscanf(stream, "%g", &tmp1)){
+            errorPrint("cscLoad: Wrong input format");
+            return EXIT_FAILURE;
+        }
+        array[i  ] = (float)(tmp1);
+        break;
+    }
+
+    return EXIT_SUCCESS;
+}
+
+/*
+ Function: csc_load
+
+ Load a csc from disk.
+
+ Fill *n*, *colptr*, *rowptr*, *values* and *dof* from *infile*.
+
+ Parameters:
+ n       - number of columns
+ colptr  - First cscd starting index of each column in *ja* and *a*
+ rowptr    - Row of each element in first CSCD
+ values  - value of each cscd in first CSCD (can be NULL)
+ dof     - Number of degrees of freedom
+ outfile - Output stream.
+
+ Return:
+ PASTIX_SUCCESS
+
+ */
+int
+csc_load( pastix_int_t  *n,
+          pastix_int_t **colptr,
+          pastix_int_t **rowptr,
+          int           *valtype,
+          void         **values,
+          int           *dof,
+          FILE          *infile )
+{
+    int  rc, ft;
+    long tmp1, tmp2;
+    pastix_int_t nnz;
+
+    /* Read the header file */
+    if (3 != fscanf(infile, "%ld %ld %d\n", &tmp1, &tmp2, &ft)) {
+        errorPrint("cscLoad:line 1: Wrong input");
+        return EXIT_FAILURE;
+    }
+
+    *n   = (pastix_int_t)tmp1;
+    *dof = (int)tmp2;
+
+    /* Read the colptr array */
+    *colptr = NULL;
+    MALLOC_INTERN(*colptr, (*n)+1, pastix_int_t);
+    assert(*colptr);
+
+    rc = readArrayOfInteger( infile, *n+1, *colptr );
+    if ( rc != EXIT_SUCCESS )
+        return rc;
+
+    /* Read the rowptr array */
+    nnz = (*colptr)[*n]-(*colptr)[0];
+    *rowptr = NULL;
+    MALLOC_INTERN(*rowptr, nnz, pastix_int_t);
+    assert(*rowptr);
+
+    rc = readArrayOfInteger( infile, nnz, *rowptr );
+    if ( rc != EXIT_SUCCESS )
+        return rc;
+
+    /* Read values if values is provided and if file contains */
+    if (values != NULL) {
+        pastix_int_t nval = nnz * (*dof) * (*dof);
+        (*values) = NULL;
+
+        switch(ft)
+        {
+        case PastixComplex64:
+            *values = malloc( nval * sizeof(pastix_complex64_t) );
+            readArrayOfComplex64( infile, nval, *values );
+            if ( rc != EXIT_SUCCESS )
+                return rc;
+            break;
+
+        case PastixComplex32:
+            *values = malloc( nval * sizeof(pastix_complex32_t) );
+            readArrayOfComplex32( infile, nval, *values );
+            if ( rc != EXIT_SUCCESS )
+                return rc;
+            break;
+
+        case PastixDouble:
+            *values = malloc( nval * sizeof(double) );
+            readArrayOfDouble( infile, nval, *values );
+            if ( rc != EXIT_SUCCESS )
+                return rc;
+            break;
+
+        case PastixFloat:
+            *values = malloc( nval * sizeof(float) );
+            readArrayOfFloat( infile, nval, *values );
+            if ( rc != EXIT_SUCCESS )
+                return rc;
+            break;
+        }
+    }
+
+    if ( valtype != NULL ) {
+        *valtype = ft;
+    }
+    return PASTIX_SUCCESS;
+}
+
+int
+cscLoad( pastix_csc_t  *csc,
+         FILE          *infile )
+{
+    int rc, dof = 1;
+
+    rc = csc_load( &(csc->n),
+                   &(csc->colptr),
+                   &(csc->rowptr),
+                   &(csc->flttype),
+                   &(csc->values),
+                   &dof,
+                   infile );
+
+    csc->gN = csc->n;
+    csc->loc2glob = NULL;
+
+    return rc;
+}
+
+
+int
+csc_save( pastix_int_t  n,
+          pastix_int_t *colptr,
+          pastix_int_t *rowptr,
+          int           ft,
+          void         *values,
+          int           dof,
+          FILE         *outfile )
+{
+    pastix_int_t i;
+
+    /* Write header N Dof FloatType */
+    fprintf( outfile, "%ld %ld %d\n",
+             (long)n, (long)dof,
+             (values == NULL) ? 0 : ft );
+
+    /* Write colptr */
+    for (i=0; i<n+1; i++)
+    {
+        fprintf(outfile, "%ld ", (long)colptr[i]);
+        if (i%4 == 3) fprintf(outfile, "\n");
+    }
+    if ((i-1)%4 !=3) fprintf(outfile, "\n");
+
+    /* Write rowptr */
+    for (i=0; i<colptr[n]-1; i++)
+    {
+        fprintf(outfile, "%ld ", (long)rowptr[i]);
+        if (i%4 == 3) fprintf(outfile, "\n");
+    }
+    if ((i-1)%4 !=3) fprintf(outfile, "\n");
+
+    /* Write the values */
+    if (values != NULL)
+    {
+        /*         for (i=0; i<(colptr[n]-1)*dof*dof; i++) */
+        /*         { */
+        /* #ifdef TYPE_COMPLEX */
+        /*             fprintf(outfile, "%lg %lg ", (double)(creal(values[i])), (double)(cimag(values[i]))); */
+        /* #else */
+        /*             fprintf(outfile, "%lg ", (double)(values[i])); */
+        /* #endif */
+        /*             if (i%4 == 3) fprintf(outfile, "\n"); */
+        /*         } */
+        /*         if ((i-1)%4 !=3) fprintf(outfile, "\n"); */
+    }
+    return PASTIX_SUCCESS;
+}
+
+int
+cscSave( pastix_csc_t *csc,
+         FILE         *outfile )
+{
+    return csc_save( csc->n,
+                     csc->colptr,
+                     csc->rowptr,
+                     csc->flttype,
+                     csc->values,
+                     1,
+                     outfile );
+}

z_csc_utils.h

+/**
+ *  PaStiX CSC management routines.
+ *
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Pierre Ramet
+ * @author Xavier Lacoste
+ * @date 2011-11-11
+ * @precisions normal z -> c d s
+ *
+ **/
+#ifndef CSC_UTILS_H
+#define CSC_UTILS_H
+
+
+/**
+ *
+ * @ingroup csc_utils
+ *
+ * Modify the CSC to a symetric graph one.
+ * Don't use it on a lower symetric CSC
+ * it would give you all the CSC upper + lower.
+ *
+ * @param[in]  n     Number of columns/vertices
+ * @param[in]  ia    Starting index of each column in *ja* and *a*
+ * @param[in]  ja    Row index of each element
+ * @param[in]  a     Value of each element,can be NULL
+ * @param[out] newn  New number of column
+ * @param[out] newia Starting index of each column in *ja* and *a*
+ * @param[out] newja Row index of each element
+ * @param[out] newa  Value of each element,can be NULL
+ **/
+int z_csc_symgraph ( pastix_int_t               n,
+                     const pastix_int_t        *ia,
+                     const pastix_int_t        *ja,
+                     const pastix_complex64_t  *a,
+                     pastix_int_t              *newn,
+                     pastix_int_t             **newia,
+                     pastix_int_t             **newja,
+                     pastix_complex64_t       **newa);
+
+
+
+
+/**
+ *
+ * @ingroup csc_utils
+ *
+ * Modify the CSC to a symetric graph one.
+ * Don't use it on a lower symetric CSC
+ * it would give you all the CSC upper + lower.
+ *
+ * Internal function.
+ *
+ * @param[in]  n           Number of columns/vertices
+ * @param[in]  ia          Starting index of each column in *ja* and *a*
+ * @param[in]  ja          Row index of each element
+ * @param[in]  a           Value of each element,can be NULL
+ * @param[out] newn        New number of column
+ * @param[out] newia       Starting index of each column in *ja* and *a*
+ * @param[out] newja       Row index of each element
+ * @param[out] newa        Value of each element,can be NULL
+ * @param[in]  malloc_flag flag to indicate if function call is intern to pastix
+ *                         or extern.
+ **/
+int z_csc_symgraph_int ( pastix_int_t               n,
+                         const pastix_int_t        *ia,
+                         const pastix_int_t        *ja,
+                         const pastix_complex64_t  *a,
+                         pastix_int_t              *newn,
+                         pastix_int_t             **newia,
+                         pastix_int_t             **newja,
+                         pastix_complex64_t       **newa,
+                         int                        malloc_flag);
+
+
+
+/**
+ *
+ * @ingroup csc_utils
+ *
+ * Supress diagonal term.
+ * After this call, *ja* can be reallocated to *ia[n] -1*.
+ *
+ * @param[in]      baseval     Initial numbering value (0 or 1).
+ * @param[in]      n           Number of columns/vertices
+ * @param[in,out]  ia          Starting index of each column in *ja* and *a*
+ * @param[in,out]  ja          Row index of each element
+ * @param[in,out]  a           Value of each element,can be NULL
+ **/
+void z_csc_noDiag( pastix_int_t        baseval,
+                   pastix_int_t        n,
+                   pastix_int_t       *ia,
+                   pastix_int_t       *ja,
+                   pastix_complex64_t *a);
+
+/**
+ *
+ * @ingroup csc_utils
+ *
+ * Check if the csc contains doubles and if correct if asked
+ *
+ * Assumes that the CSC is sorted.
+ *
+ * Assumes that the CSC is Fortran numeroted (base 1)
+ *
+ * @param[in]     n          Size of the matrix.
+ * @param[in,out] colptr     Index in *rows* and *values* of the first element
+ *                           of each column
+ * @param[in,out] rows       row of each element
+ * @param[in,out] values     value of each element
+ * @param[in]     dof        Number of degrees of freedom
+ * @param[in]     flag       Indicate if user wants correction (<API_BOOLEAN>)
+ * @param[in]     flagalloc  indicate if allocation on CSC uses internal malloc.
+ *
+ *
+ * @Returns:
+ *   API_YES - If the matrix contained no double or was successfully corrected.
+ *   API_NO  - Otherwise.
+*/
+
+int z_csc_check_doubles(pastix_int_t         n,
+                        pastix_int_t        *colptr,
+                        pastix_int_t       **rows,
+                        pastix_complex64_t **values,
+                        int                  dof,
+                        int                  flag,
+                        int                  flagalloc);
+
+
+/**
+ *
+ * @ingroup csc_utils
+ *
+ * Check if the CSC graph is symetric.
+ *
+ *   For all local column C,
+ *
+ *   For all row R in the column C,
+ *
+ *   We look in column R if we have the row number C.
+ *
+ *   If we can correct we had missing non zeros.
+ *
+ *   Assumes that the CSC is Fortran numbered (1 based).
+ *
+ *   Assumes that the matrix is sorted.
+ *
+ * @param[in]     n         Number of local columns
+ * @param[in,out] colptr    Starting index of each columns in *ja*
+ * @param[in,out] rows      Row of each element.
+ * @param[in,out] values    Value of each element.
+ * @param[in]     correct   Flag indicating if we can correct the symmetry.
+ * @param[in]     alloc     indicate if allocation on CSC uses internal malloc.
+ * @param[in]     dof       Number of degrees of freedom.
+*/
+int z_csc_checksym(pastix_int_t         n,
+                   pastix_int_t        *colptr,
+                   pastix_int_t       **rows,
+                   pastix_complex64_t **values,
+                   int                  correct,
+                   int                  alloc,
+                   int                  dof);
+
+void z_csc_colPerm(pastix_int_t n, pastix_int_t *ia, pastix_int_t *ja, pastix_complex64_t *a, pastix_int_t *cperm);
+void z_csc_colScale(pastix_int_t n, pastix_int_t *ia, pastix_int_t *ja, pastix_complex64_t *a, pastix_complex64_t *dcol);
+void z_csc_rowScale(pastix_int_t n, pastix_int_t *ia, pastix_int_t *ja, pastix_complex64_t *a, pastix_complex64_t *drow);
+
+void z_csc_sort(pastix_int_t n, pastix_int_t *ia, pastix_int_t *ja, pastix_complex64_t *a, pastix_int_t ndof);
+void z_csc_Fnum2Cnum(pastix_int_t *ja, pastix_int_t *ia, pastix_int_t n);
+void z_csc_Cnum2Fnum(pastix_int_t *ja, pastix_int_t *ia, pastix_int_t n);
+
+/*
+  Function: CSC_buildZerosAndNonZerosGraphs
+
+  Separate a graph in two graphs, following
+  wether the diagonal term of a column is null or not.
+
+  Parameters:
+    n, colptr, rows, values  - The initial CSC
+    n_nz, colptr_nz, rows_nz - The graph of the non-null diagonal part.
+    n_z, colptr_z, rows_z    - The graph of the null diagonal part.
+    perm                     - Permutation to go from the first graph to
+                               the one composed of the two graph concatenated.
+    revperm                  - Reverse permutation tabular.
+    criteria                 - Value beside which a number is said null.
+*/
+int z_csc_buildZerosAndNonZerosGraphs(pastix_int_t     n,
+                                    pastix_int_t    *colptr,
+                                    pastix_int_t    *rows,
+                                    pastix_complex64_t  *values,
+                                    pastix_int_t    *n_nz,
+                                    pastix_int_t   **colptr_nz,
+                                    pastix_int_t   **rows_nz,
+                                    pastix_int_t    *n_z,
+                                    pastix_int_t   **colptr_z,
+                                    pastix_int_t   **rows_z,
+                                    pastix_int_t    *perm,
+                                    pastix_int_t    *revperm,
+                                    double  criteria);
+
+/*
+  Function: CSC_isolate
+
+  Isolate a list of unknowns at the end of the CSC.
+
+  Parameters:
+    n            - Number of columns.
+    colptr       - Index of first element of each column in *ia*.
+    rows         - Rows of each non zeros.
+    n_isolate    - Number of unknow to isolate.
+    isolate_list - List of unknown to isolate.
+*/
+int z_csc_isolate(pastix_int_t     n,
+                pastix_int_t    *colptr,
+                pastix_int_t    *rows,
+                pastix_int_t     n_isolate,
+                pastix_int_t    *isolate_list,
+                pastix_int_t    *perm,
+                pastix_int_t    *revperm);
+
+
+/*
+  Function: csc_save
+
+  Save a csc on disk.
+
+  Parameters:
+    n       - number of columns
+    colptr  - First cscd starting index of each column in *ja* and *a*
+    rows    - Row of each element in first CSCD
+    values  - value of each cscd in first CSCD (can be NULL)
+    dof     - Number of degrees of freedom
+    outfile - Output stream.
+
+  Return:
+    NO_ERR
+
+*/
+int z_csc_save(pastix_int_t      n,
+             pastix_int_t    * colptr,
+             pastix_int_t    * rows,
+             pastix_complex64_t  * values,
+             int      dof,
+             FILE   * outfile);
+/*
+  Function: csc_load
+
+  Load a csc from disk.
+
+  Fill *n*, *colptr*, *rows*, *values* and *dof* from *infile*.
+
+  Parameters:
+    n       - number of columns
+    colptr  - First cscd starting index of each column in *ja* and *a*
+    rows    - Row of each element in first CSCD
+    values  - value of each cscd in first CSCD (can be NULL)
+    dof     - Number of degrees of freedom
+    outfile - Output stream.
+
+  Return:
+    NO_ERR
+
+*/
+int z_csc_load(pastix_int_t    *  n,
+             pastix_int_t    ** colptr,
+             pastix_int_t    ** rows,
+             pastix_complex64_t  ** values,
+             int    *  dof,
+             FILE   *  infile);
+
+#endif /* CSC_UTILS_H */

z_cscd_utils.h

+/**
+ *  PaStiX distributed CSC management routines.
+ *
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Pierre Ramet
+ * @author Xavier Lacoste
+ * @date 2011-11-11
+ * @precisions normal z -> c d s
+ *
+ **/
+#ifndef Z_CSCD_UTILS_H
+#define Z_CSCD_UTILS_H
+
+
+/*
+ * Enum: CSCD_OPERATIONS
+ *
+ * Operation when adding CSCD
+ *
+ * CSCD_ADD  - Add coefficient values.
+ * CSCD_KEEP - Keep value from first CSCD.
+ * CSCD_MAX  - Keep maximum of first and second CSCD.
+ * CSCD_MIN  - Keep minimum of first and second CSCD.
+ * CSCD_OVW  - Overwrite with second CSCD value.
+ */
+enum CSCD_OPERATIONS {
+  CSCD_ADD,
+  CSCD_KEEP,
+  CSCD_MAX,
+  CSCD_MIN,
+  CSCD_OVW
+};
+typedef enum CSCD_OPERATIONS CSCD_OPERATIONS_t;
+
+/*
+ * Enum: CSC_DISPATCH_OP
+ *
+ * Operation when dispatching the csc into a cscd
+ *
+ * CSC_DISP_SIMPLE - Reparts linearly the columns over the proc.
+ * CSC_DISP_CYCLIC - Reparts cyclicly the columns over the proc.
+ *
+ */
+enum CSC_DISPATCH_OP {
+  CSC_DISP_SIMPLE,
+  CSC_DISP_CYCLIC
+};
+typedef enum CSC_DISPATCH_OP CSCD_DISPATCH_OP_t;
+
+/* Section: Functions */
+/*
+ *  Function: z_csc_dispatch
+ *
+ *  Distribute a CSC to a CSCD
+ *
+ *  Parameters:
+ *     gN                - global number of columns
+ *     gcolptr           - global starting index of each column in grows ans gavals.
+ *     grows             - global rows of each element.
+ *     gavals            - global values of each element.
+ *     gperm             - global permutation tabular.
+ *     ginvp             - global reverse permutation tabular.
+ *     lN                - local number of columns (output).
+ *     lcolptr           - starting index of each local column (output).
+ *     lrowptr           - row number of each local element (output).
+ *     lavals            - values of each local element (output).
+ *     lrhs              - local part of the right hand side (output).
+ *     lperm             - local part of the permutation tabular (output).
+ *     loc2glob          - global numbers of local columns (before permutation).
+ *     dispatch          - choose how to dispatch the csc
+ *     pastix_comm       - PaStiX MPI communicator.
+ */
+void z_csc_dispatch(pastix_int_t  gN, pastix_int_t *  gcolptr, pastix_int_t *  grow, pastix_complex64_t *  gavals,
+                  pastix_complex64_t *  grhs, pastix_int_t *  gperm, pastix_int_t *  ginvp,
+                  pastix_int_t *lN, pastix_int_t ** lcolptr, pastix_int_t ** lrow, pastix_complex64_t ** lavals,
+                  pastix_complex64_t ** lrhs, pastix_int_t ** lperm,
+                  pastix_int_t **loc2glob, int dispatch, MPI_Comm pastix_comm);
+
+/*
+ * Function: csc_cyclic_distribution
+ *
+ * Distribute the CSC cyclicaly.
+ *
+ * Parameters:
+ *   column      - column number to distribute
+ *   columnnbr   - Number of colmuns.
+ *   pastix_comm - PaStiX MPI communicator
+ *
+ * Return:
+ *   owner of the column (column%commSize)
+ */
+pastix_int_t z_csc_cyclic_distribution(pastix_int_t column, pastix_int_t columnnbr, MPI_Comm pastix_comm);
+
+/*
+ * Function: csc_simple_distribution
+ *
+ * Distribute the CSC.
+ * First columns are for first proc and so on.
+ *
+ * Parameters:
+ *   column      - column number to distribute
+ *   columnnbr   - Number of colmuns.
+ *   pastix_comm - PaStiX MPI communicator
+ *
+ * Return:
+ *   owner of the column (column/commSize)
+ */
+pastix_int_t z_csc_simple_distribution(pastix_int_t column, pastix_int_t columnnbr, MPI_Comm pastix_comm);
+
+/*
+ * Function: cscd_symgraph
+ *
+ * Check if the CSCD graph is symetric.
+ *
+ * Parameters:
+ *   n           - Number of local columns
+ *   ia          - Starting index of each columns in *ja* and *a*
+ *   ja          - Row of each element.
+ *   a           - Values of each element.
+ *   newn        - New number of local columns
+ *   newia       - Starting index of each columns in *newja* and *newa*
+ *   newja       - Row of each element.
+ *   newa        - Values of each element.
+ *   l2g         - global number of each local column.
+ *   malloc_flag - flag to indicate if function call is intern to pastix or extern.
+ */
+int z_cscd_symgraph(pastix_int_t      n, const pastix_int_t *ia, const pastix_int_t *ja, const pastix_complex64_t *a,
+                    pastix_int_t * newn, pastix_int_t ** newia, pastix_int_t ** newja, pastix_complex64_t ** newa,
+                    pastix_int_t *     l2g,  MPI_Comm comm);
+/*
+ * Function: cscd_addlocal
+ *
+ * Add second cscd to first cscd into third cscd (unallocated)
+ *
+ * Parameters:
+ *   n           - First cscd size
+ *   ia          - First cscd starting index of each column in *ja* and *a*
+ *   ja          - Row of each element in first CSCD
+ *   a           - value of each cscd in first CSCD (can be NULL)
+ *   l2g         - local 2 global column numbers for first cscd
+ *   addn        - CSCD to add size
+ *   addia       - CSCD to add starting index of each column in *addja* and *adda*
+ *   addja       - Row of each element in second CSCD
+ *   adda        - value of each cscd in second CSCD (can be NULL -> add 0)
+ *   addl2g      - local 2 global column numbers for second cscd
+ *   newn        - new cscd size (same as first)
+ *   newia       - CSCD to add starting index of each column in *newja* and *newwa*
+ *   newja       - Row of each element in third CSCD
+ *   newa        - value of each cscd in third CSCD
+ *   OP          - Operation to manage common CSCD coefficients.
+ *   dof         - Number of degrees of freedom.
+ */
+
+int z_cscd_addlocal(pastix_int_t   n   , pastix_int_t *  ia   , pastix_int_t *  ja   , pastix_complex64_t *  a   , pastix_int_t * l2g,
+      pastix_int_t   addn, pastix_int_t *  addia, pastix_int_t *  addja, pastix_complex64_t *  adda, pastix_int_t * addl2g,
+      pastix_int_t * newn, pastix_int_t ** newia, pastix_int_t ** newja, pastix_complex64_t ** newa, CSCD_OPERATIONS_t OP, int dof);
+
+/**
+ *   Function: csc2cscd
+ *
+ *   Transform a csc to a cscd.
+ *   Allocate the CSCD.
+ *   If grhs == NULL forget right hand side part.
+ *   If gperm == NULL forget permutation and reverse permutation part.
+ *
+ *   Parameters:
+ *     gN       - global number of columns
+ *     gcolptr  - global starting index of each column in grows ans gavals.
+ *     grows    - global rows of each element.
+ *     gavals   - global values of each element.
+ *     gperm    - global permutation tabular.
+ *     ginvp    - global reverse permutation tabular.
+ *     lN       - local number of columns.
+ *     lcolptr  - starting index of each local column.
+ *     lrowptr  - row number of each local element.
+ *     lavals   - values of each local element.
+ *     lrhs     - local part of the right hand side (output).
+ *     lperm    - local part of the permutation tabular (output).
+ *     linvp    - local part of the reverse permutation tabular (output).
+ *     loc2glob - global numbers of local columns (before permutation).
+ */
+void  z_csc2cscd(pastix_int_t gN, pastix_int_t *  gcolptr, pastix_int_t *  grow, pastix_complex64_t *  gavals,
+               pastix_complex64_t *  grhs, pastix_int_t *  gperm, pastix_int_t *  ginvp,
+               pastix_int_t lN, pastix_int_t ** lcolptr, pastix_int_t ** lrow, pastix_complex64_t ** lavals,
+               pastix_complex64_t ** lrhs, pastix_int_t ** lperm, pastix_int_t ** linvp,
+               pastix_int_t *loc2glob);
+
+/**
+ *   Function: cscd2csc
+ *
+ *   Transform a cscd to a csc.
+ *   colptr2, row2, avals2, rhs2, perm2, invp2 are allocated here.
+ *
+ *   Parameters:
+ *      lN          - number of local column.
+ *      lcolptr     - starting index of each local column in row and avals.
+ *      lrow        _ row number of each local element.
+ *      lavals      - values of each local element.
+ *      lrhs        - local part of the right hand side.
+ *      lperm       - local part of the permutation tabular.
+ *      linvp       - local part of the reverse permutation tabular.
+ *      gN          - global number of columns (output).
+ *      gcolptr     - starting index of each column in row2 and avals2 (output).
+ *      grow        - row number of each element (output).
+ *      gavals      - values of each element (output).
+ *      grhs        - global right hand side (output).
+ *      gperm       - global permutation tabular (output).
+ *      ginvp       - global reverse permutation tabular (output).
+ *      loc2glob    - global number of each local column.
+ *      pastix_comm - PaStiX MPI communicator.
+ *      ndof        - Number of degree f freedom by node.
+ *
+ */
+
+void  z_cscd2csc(pastix_int_t  lN, pastix_int_t *  lcolptr, pastix_int_t * lrow, pastix_complex64_t * lavals,
+               pastix_complex64_t * lrhs, pastix_int_t * lperm, pastix_int_t * linvp,
+               pastix_int_t *gN, pastix_int_t ** gcolptr, pastix_int_t **grow, pastix_complex64_t **gavals,
+               pastix_complex64_t **grhs, pastix_int_t **gperm, pastix_int_t **ginvp,
+               pastix_int_t *loc2glob, MPI_Comm pastix_comm, pastix_int_t ndof);
+
+/*
+ * Function: cscd_redispatch
+ *
+ * Redistribute the first cscd into a new one using *dl2g*.
+ *
+ * - gather all new loc2globs on all processors.
+ * - allocate *dia*, *dja* and *da*.
+ * - Create new CSC for each processor and send it.
+ * - Merge all new CSC to the new local CSC with <cscd_addlocal_int>.
+ *
+ * If communicator size is one, check that n = dn and
+ * l2g = dl2g and simply create a copy of the first cscd.
+ *
+ * Parameters:
+ *   n           - Number of local columns
+ *   ia          - First cscd starting index of each column in *ja* and *a*
+ *   ja          - Row of each element in first CSCD
+ *   a           - value of each cscd in first CSCD (can be NULL)
+ *   rhs         - right-hand-side member corresponding to the first CSCD
+ *                 (can be NULL)
+ *   nrhs        - number of right-hand-side.
+ *   l2g         - local 2 global column numbers for first cscd
+ *   dn          - Number of local columns
+ *   dia         - New cscd starting index of each column in *ja* and *a*
+ *   dja         - Row of each element in new CSCD
+ *   da          - value of each cscd in new CSCD
+ *   rhs         - right-hand-side member corresponding to the new CSCD
+ *   dl2g        - local 2 global column numbers for new cscd
+ *   comm        - MPI communicator
+ *
+ * Returns:
+ *   EXIT_SUCCESS - If all goes well
+ *   EXIT_FAILURE - If commsize = 1 and *n* != *dn* or *l2g* != *dl2g*.
+ */
+int z_cscd_redispatch(pastix_int_t   n, pastix_int_t *   ia, pastix_int_t *   ja, pastix_complex64_t *   a,
+                    pastix_complex64_t *  rhs,  pastix_int_t nrhs, pastix_int_t *   l2g,
+                    pastix_int_t  dn, pastix_int_t ** dia, pastix_int_t ** dja, pastix_complex64_t ** da,
+                    pastix_complex64_t ** drhs,  pastix_int_t *  dl2g,
+                    MPI_Comm comm, pastix_int_t dof);
+
+/*
+ *  Function: cscd_save
+ *
+ *  save a distributed csc to disk.
+ *  files are called $(filename) and $(filename)$(RANK)
+ *  if filename is NULL then filename = cscd_matrix.
+ *
+ *  file filename contains the number of processors/files
+ *  on first line. Then each line contain the name of each file
+ *  (here $(filename)$(RANK)).
+ *
+ *
+ *
+ *  Parameters:
+ *    n           - Number of local columns
+ *    ia          - First cscd starting index of each column in *ja* and *a*
+ *    ja          - Row of each element in first CSCD
+ *    a           - value of each cscd in first CSCD (can be NULL)
+ *    rhs         - Right hand side.
+ *    l2g         - local 2 global column numbers for first cscd
+ *    dof         - Number of degrees of freedom
+ *    filename    - name of the files.
+ *    comm        - MPI communicator
+ */
+int z_cscd_save(pastix_int_t n, pastix_int_t *ia, pastix_int_t *ja, pastix_complex64_t * a, pastix_complex64_t * rhs, pastix_int_t* l2g,
+              int dof, const char * filename, MPI_Comm comm);
+
+/*
+ *  Function: cscd_load
+ *
+ *  Loads a distributed csc from disk.
+ *  if filename is NULL then filename = cscd_matrix.
+ *
+ *  Parameters:
+ *    n           - Number of local columns
+ *    ia          - First cscd starting index of each column in *ja* and *a*
+ *    ja          - Row of each element in first CSCD
+ *    a           - value of each cscd in first CSCD (can be NULL)
+ *    rhs         - Right hand side.
+ *    l2g         - local 2 global column numbers for first cscd
+ *    filename    - name of the files.
+ *    comm        - MPI communicator
+ */
+int z_cscd_load(pastix_int_t *n, pastix_int_t ** ia, pastix_int_t ** ja, pastix_complex64_t ** a, pastix_complex64_t ** rhs, pastix_int_t ** l2g,
+              const char * filename, MPI_Comm mpi_comm);
+
+#endif /* Z_CSCD_UTILS_H */

z_cscd_utils_fortran.c

+/**
+ *  PaStiX CSC management routines.
+ *
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Pierre Ramet
+ * @author Xavier Lacoste
+ * @date 2011-11-11
+ * @precisions normal z -> c d s
+ *
+ **/
+/*
+  File: pastix_fortran.c
+
+  Interface to the PaStiX API functions.
+
+ */
+#include "common.h"
+#include "z_cscd_utils.h"
+
+
+#define FORTRAN_NAME(nu,nl,pl,pc)   \
+  void nl pl;                       \
+  void nu pl                        \
+  { nl pc; }                        \
+  void nl ## _ pl                   \
+  { nl pc; }                        \
+  void nl ## __ pl                  \
+  { nl pc; }
+/*
+  Struct: csc_data_
+
+  Contains the new CSCD
+
+*/
+struct csc_data_ {
+  pastix_int_t        n;
+  pastix_int_t       *colptr;
+  pastix_int_t       *rows;
+  pastix_complex64_t *values;
+  pastix_complex64_t *rhs;
+  pastix_int_t        nrhs;
+  pastix_int_t       *perm;
+  pastix_int_t       *l2g;
+  pastix_int_t        dof;
+};
+
+/*
+  Typedef: csc_data_t
+
+  Type coresponding to the struct <csc_data_>
+*/
+typedef struct csc_data_ csc_data_t;
+
+FORTRAN_NAME(Z_CSC_DISPATCH_FORTRAN,
+             z_csc_dispatch_fortran,
+             (csc_data_t         **csc_data,
+              pastix_int_t        *gN,
+              pastix_int_t        *gcolptr,
+              pastix_int_t        *grow,
+              pastix_complex64_t  *gavals,
+              pastix_complex64_t  *grhs,
+              pastix_int_t        *gperm,
+              pastix_int_t        *ginvp,
+              int                 *dispatch,
+              pastix_int_t        *newn,
+              pastix_int_t        *newnnz,
+              MPI_Fint            *fortran_comm),
+             (csc_data,
+              gN,
+              gcolptr,
+              grow,
+              gavals,
+              grhs,
+              gperm,
+              ginvp,
+              dispatch,
+              newn,
+              newnnz,
+              fortran_comm))
+
+void z_csc_dispatch_fortran(csc_data_t         **csc_data,
+                            pastix_int_t        *gN,
+                            pastix_int_t        *gcolptr,
+                            pastix_int_t        *grow,
+                            pastix_complex64_t  *gavals,
+                            pastix_complex64_t  *grhs,
+                            pastix_int_t        *gperm,
+                            pastix_int_t        *ginvp,
+                            int                 *dispatch,
+                            pastix_int_t        *newn,
+                            pastix_int_t        *newnnz,
+                            MPI_Fint            *fortran_comm)
+{
+  MPI_Comm        pastix_comm;
+  (void)fortran_comm;
+
+  pastix_comm = MPI_Comm_f2c(*fortran_comm);
+  MALLOC_INTERN(*csc_data, 1, struct csc_data_);
+  (*csc_data)->n      = 0;
+  (*csc_data)->colptr = NULL;
+  (*csc_data)->rows   = NULL;
+  (*csc_data)->values = NULL;
+  (*csc_data)->rhs    = NULL;
+  (*csc_data)->perm   = NULL;
+  (*csc_data)->l2g    = NULL;
+
+  z_csc_dispatch(*gN, gcolptr, grow, gavals, grhs, gperm, ginvp,
+                 &((*csc_data)->n), &((*csc_data)->colptr), &((*csc_data)->rows), &((*csc_data)->values),
+               &((*csc_data)->rhs), &((*csc_data)->perm),
+               &((*csc_data)->l2g), *dispatch, pastix_comm);
+
+  *newn   = (*csc_data)->n;
+  *newnnz = (*csc_data)->colptr[(*csc_data)->n]-1;
+
+}
+
+
+
+FORTRAN_NAME(Z_CSC_DISPATCH_FORTRAN_END,
+             z_csc_dispatch_fortran_end,
+             (csc_data_t         **csc_data,
+              pastix_int_t        *lcolptr,
+              pastix_int_t        *lrow,
+              pastix_complex64_t  *lavals,
+              pastix_complex64_t  *lrhs,
+              pastix_int_t        *lperm,
+              pastix_int_t        *l2g),
+             (csc_data,
+              lcolptr,
+              lrow,
+              lavals,
+              lrhs,
+              lperm,
+              l2g))
+void z_csc_dispatch_fortran_end(csc_data_t         **csc_data,
+                                pastix_int_t        *lcolptr,
+                                pastix_int_t        *lrow,
+                                pastix_complex64_t  *lavals,
+                                pastix_complex64_t  *lrhs,
+                                pastix_int_t       *lperm,
+                                pastix_int_t        *l2g)
+{
+  pastix_int_t nnz = 0;
+  if ((*csc_data)->colptr != NULL) {
+    nnz = (*csc_data)->colptr[(*csc_data)->n]-1;
+    memcpy(lcolptr, (*csc_data)->colptr, (1+(*csc_data)->n)*sizeof(pastix_int_t));
+    free((*csc_data)->colptr);
+  }
+  if ((*csc_data)->rows != NULL) {
+    memcpy(lrow,    (*csc_data)->rows,   nnz*sizeof(pastix_int_t));
+    free((*csc_data)->rows);
+  }
+  if ((*csc_data)->values != NULL) {
+    memcpy(lavals,  (*csc_data)->values,   nnz*sizeof(pastix_complex64_t));
+    free((*csc_data)->values);
+  }
+  if ((*csc_data)->rhs != NULL) {
+    memcpy(lrhs, (*csc_data)->rhs, (*csc_data)->n*sizeof(pastix_complex64_t));
+    free((*csc_data)->rhs);
+  }
+  if ((*csc_data)->perm != NULL) {
+    memcpy(lperm, (*csc_data)->perm, (*csc_data)->n*sizeof(pastix_int_t));
+    free((*csc_data)->perm);
+  }
+  if ((*csc_data)->l2g != NULL) {
+    memcpy(l2g, (*csc_data)->l2g, (*csc_data)->n*sizeof(pastix_int_t));
+    free((*csc_data)->l2g);
+  }
+}
+
+FORTRAN_NAME(Z_CSCD_REDISPATCH_FORTRAN,
+             z_cscd_redispatch_fortran,
+             (csc_data_t         **csc_data,
+              pastix_int_t        *n,
+              pastix_int_t        *ia,
+              pastix_int_t        *ja,
+              pastix_complex64_t  *a,
+              pastix_complex64_t  *rhs,
+              pastix_int_t        *nrhs,
+              pastix_int_t        *l2g,
+              pastix_int_t        *newn,
+              pastix_int_t        *newl2g,
+              pastix_int_t        *newnnz,
+              MPI_Fint            *fortran_comm,
+              pastix_int_t        *dof,
+              int                 *ierr),
+             (csc_data,
+              n,
+              ia,
+              ja,
+              a,
+              rhs,
+              nrhs,
+              l2g,
+              newn,
+              newl2g,
+              newnnz,
+              fortran_comm,
+              dof,
+              ierr))
+void z_cscd_redispatch_fortran(csc_data_t         **csc_data,
+                               pastix_int_t        *n,
+                               pastix_int_t        *ia,
+                               pastix_int_t        *ja,
+                               pastix_complex64_t  *a,
+                               pastix_complex64_t  *rhs,
+                               pastix_int_t        *nrhs,
+                               pastix_int_t        *l2g,
+                               pastix_int_t        *newn,
+                               pastix_int_t        *newl2g,
+                               pastix_int_t        *newnnz,
+                               MPI_Fint            *fortran_comm,
+                               pastix_int_t        *dof,
+                               int                 *ierr)
+{
+
+  MPI_Comm        pastix_comm;
+  (void)newnnz; (void)fortran_comm;
+
+  pastix_comm = MPI_Comm_f2c(*fortran_comm);
+  MALLOC_INTERN(*csc_data, 1, struct csc_data_);
+  (*csc_data)->n      = 0;
+  (*csc_data)->colptr = NULL;
+  (*csc_data)->rows   = NULL;
+  (*csc_data)->values = NULL;
+  (*csc_data)->rhs    = NULL;
+  (*csc_data)->nrhs   = *nrhs;
+  (*csc_data)->perm   = NULL;
+  (*csc_data)->l2g    = NULL;
+  (*csc_data)->dof    = *dof;
+
+  *ierr = z_cscd_redispatch(*n,    ia, ja, a, rhs,  *nrhs, l2g,
+                            *newn, &((*csc_data)->colptr), &((*csc_data)->rows), &((*csc_data)->values),  &((*csc_data)->rhs), newl2g,
+                          pastix_comm, (*csc_data)->dof);
+}
+
+FORTRAN_NAME(Z_CSCD_REDISPATCH_FORTRAN_END,
+             z_cscd_redispatch_fortran_end,
+             (csc_data_t         **csc_data,
+              pastix_int_t        *dcolptr,
+              pastix_int_t        *drow,
+              pastix_complex64_t  *davals,
+              pastix_complex64_t  *drhs),
+             (csc_data,
+              dcolptr,
+              drow,
+              davals,
+              drhs))
+
+
+void z_cscd_redispatch_fortran_end(csc_data_t         **csc_data,
+                                   pastix_int_t        *dcolptr,
+                                   pastix_int_t        *drow,
+                                   pastix_complex64_t  *davals,
+                                   pastix_complex64_t  *drhs)    
+{
+  pastix_int_t nnz = 0;
+  if ((*csc_data)->colptr != NULL) {
+    nnz = (*csc_data)->colptr[(*csc_data)->n]-1;
+    memcpy(dcolptr, (*csc_data)->colptr, (1+(*csc_data)->n)*sizeof(pastix_int_t));
+    free((*csc_data)->colptr);
+  }
+  if ((*csc_data)->rows != NULL) {
+    memcpy(drow,    (*csc_data)->rows,   nnz*sizeof(pastix_int_t));
+    free((*csc_data)->rows);
+  }
+  if ((*csc_data)->values != NULL) {
+    memcpy(davals,  (*csc_data)->values,
+           (*csc_data)->dof*(*csc_data)->dof*nnz*sizeof(pastix_complex64_t));
+    free((*csc_data)->values);
+  }
+  if ((*csc_data)->rhs != NULL) {
+    memcpy(drhs, (*csc_data)->rhs, (*csc_data)->nrhs*(*csc_data)->n*sizeof(pastix_complex64_t));
+    free((*csc_data)->rhs);
+  }
+  memFree_null(*csc_data);
+}

z_cscd_utils_intern.h

+/**
+ *  PaStiX distributed CSC management routines.
+ *
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Pierre Ramet
+ * @author Xavier Lacoste
+ * @date 2011-11-11
+ * @precisions normal z -> c d s
+ *
+ **/
+#include "z_cscd_utils.h"
+#ifndef CSCD_UTILS_INTERN_H
+#define CSCD_UTILS_INTERN_H
+int z_cscd_addlocal_int(pastix_int_t   n   , const pastix_int_t *  ia   , const pastix_int_t *  ja   , const pastix_complex64_t *  a   , const pastix_int_t * l2g,
+                        pastix_int_t   addn, const pastix_int_t *  addia,       pastix_int_t *  addja, const pastix_complex64_t *  adda, const pastix_int_t * addl2g,
+                        pastix_int_t * newn, pastix_int_t ** newia, pastix_int_t ** newja, pastix_complex64_t ** newa,
+                        CSCD_OPERATIONS_t OP, int dof, int malloc_flag);
+
+/*
+ * Function: cscd_redispatch_int
+ *
+ * Redistribute the first cscd into a new one using *dl2g*.
+ *
+ * - gather all new loc2globs on all processors.
+ * - allocate *dia*, *dja* and *da*.
+ * - Create new CSC for each processor and send it.
+ * - Merge all new CSC to the new local CSC with <cscd_addlocal_int>.
+ *
+ * If communicator size is one, check that n = dn and
+ * l2g = dl2g and simply create a copy of the first cscd.
+ *
+ * Parameters:
+ *   n           - Number of local columns
+ *   ia          - First cscd starting index of each column in *ja* and *a*
+ *   ja          - Row of each element in first CSCD
+ *   a           - value of each cscd in first CSCD (can be NULL)
+ *   rhs         - right-hand-side member corresponding to the first CSCD (can be NULL)
+ *   nrhs        - number of right-hand-side.
+ *   l2g         - local 2 global column numbers for first cscd
+ *   dn          - Number of local columns
+ *   dia         - New cscd starting index of each column in *ja* and *a*
+ *   dja         - Row of each element in new CSCD
+ *   da          - value of each cscd in new CSCD
+ *   rhs         - right-hand-side member corresponding to the new CSCD
+ *   dl2g        - local 2 global column numbers for new cscd
+ *   malloc_flag - Internal (API_YES) or external (API_NO) malloc use.
+ *   comm        - MPI communicator
+ *
+ * Returns:
+ *   EXIT_SUCCESS - If all goes well
+ *   EXIT_FAILURE - If commsize = 1 and *n* != *dn* or *l2g* != *dl2g*.
+ */
+int z_cscd_redispatch_int(pastix_int_t   n, pastix_int_t *   ia, pastix_int_t *   ja, pastix_complex64_t *   a, pastix_complex64_t *  rhs,  pastix_int_t nrhs, pastix_int_t *   l2g,
+                        pastix_int_t  dn, pastix_int_t ** dia, pastix_int_t ** dja, pastix_complex64_t ** da, pastix_complex64_t ** drhs, pastix_int_t *  dl2g,
+                        int  malloc_flag, MPI_Comm comm, pastix_int_t dof);
+
+int z_cscd_symgraph_int(pastix_int_t      n, const pastix_int_t *ia, const pastix_int_t *ja, const pastix_complex64_t *a,
+                        pastix_int_t * newn, pastix_int_t **  newia, pastix_int_t **  newja, pastix_complex64_t ** newa,
+                        pastix_int_t *  l2g, MPI_Comm comm, int malloc_flag);
+
+
+/*
+  Function: cscd_build_g2l
+
+  Construct global to local tabular containing local number of global columns
+  if one column is local, and -owner if column is not local.
+
+  For i in 0, gN
+     g2l[i] = i local number if i is local
+     g2l[i] = -p if p is the owner of the column i
+
+  Parameters:
+    n        - Number of local columns
+    colptr   - Starting index of each columns in *ja*
+    rows     - Row of each element.
+    values   - Value of each element.
+    l2g      - global number of each local column.
+    correct  - Flag indicating if we can correct the symmetry.
+    dof      - Number of degrees of freedom.
+    comm     - MPI communicator
+ */
+int z_cscd_build_g2l( pastix_int_t   ncol,
+                      pastix_int_t  *loc2glob,
+                      MPI_Comm       comm,
+                      pastix_int_t  *gN,
+                      pastix_int_t **g2l);
+/*
+   Function: cscd_noDiag
+
+   Removes diagonal elements from a CSCD.
+   *ja* and *a* can be reallocated to
+   ia[n]-1 elements after this call.
+
+   Parameters:
+     n           - Number of local columns
+     ia          - First cscd starting index of each column in *ja* and *a*
+     ja          - Row of each element in first CSCD
+     a           - value of each cscd in first CSCD (can be NULL)
+     l2g         - local 2 global column numbers for first cscd
+
+*/
+int z_cscd_noDiag(pastix_int_t n, pastix_int_t *ia, pastix_int_t *ja, pastix_complex64_t * a, const pastix_int_t * l2g);
+
+
+/*
+  Function: cscd_checksym
+
+  Check if the CSCD graph is symetric.
+
+  Parameters:
+    n   - Number of local columns
+    ia  - Starting index of each columns in *ja*
+    ja  - Row of each element.
+    l2g - global number of each local column.
+    correct  - Flag indicating if we can correct the symmetry.
+    alloc    - indicate if allocation on CSC uses internal malloc.
+    dof      - Number of degrees of freedom.
+    comm     - MPI communicator
+*/
+int z_cscd_checksym(pastix_int_t      n,
+                  pastix_int_t     *colptr,
+                  pastix_int_t    **rows,
+                  pastix_complex64_t  **values,
+                  pastix_int_t     *l2g,
+                  int      correct,
+                  int      alloc,
+                  int      dof,
+                  MPI_Comm comm);
+
+
+/**
+ *   Function: cscd2csc_int
+ *
+ *   Transform a cscd to a csc.
+ *   colptr2, row2, avals2, rhs2, perm2, invp2 are allocated here.
+ *
+ *   External function, allocation are not of the internal type.
+ *
+ *   Parameters:
+ *      lN          - number of local column.
+ *      lcolptr     - starting index of each local column in row and avals.
+ *      lrow        _ row number of each local element.
+ *      lavals      - values of each local element.
+ *      lrhs        - local part of the right hand side.
+ *      lperm       - local part of the permutation tabular.
+ *      linvp       - Means nothing, to suppress.
+ *      gN          - global number of columns (output).
+ *      gcolptr     - starting index of each column in row2 and avals2 (output).
+ *      grow        - row number of each element (output).
+ *      gavals      - values of each element (output).
+ *      grhs        - global right hand side (output).
+ *      gperm       - global permutation tabular (output).
+ *      ginvp       - global reverse permutation tabular (output).
+ *      loc2glob    - global number of each local column.
+ *      pastix_comm - PaStiX MPI communicator.
+ *      ndof        - Number of degree of freedom per node.
+ *      intern_flag - Decide if malloc will use internal or external macros.
+ */
+void  z_cscd2csc_int(pastix_int_t  lN, pastix_int_t *  lcolptr, pastix_int_t * lrow, pastix_complex64_t * lavals,
+                   pastix_complex64_t * lrhs, pastix_int_t * lperm, pastix_int_t * linvp,
+                   pastix_int_t *gN, pastix_int_t ** gcolptr, pastix_int_t **grow, pastix_complex64_t **gavals,
+                   pastix_complex64_t **grhs, pastix_int_t **gperm, pastix_int_t **ginvp,
+                   pastix_int_t *loc2glob, MPI_Comm pastix_comm, pastix_int_t ndof, int intern_flag);
+
+/*
+  Function: cscd_redispatch_scotch
+
+  Renumber the columns to have first columns on first proc for Scotch
+
+  Parameters:
+    n           - Number of local columns
+    ia          - First cscd starting index of each column in *ja* and *a*
+    ja          - Row of each element in first CSCD
+    a           - value of each cscd in first CSCD (can be NULL)
+    l2g         - local 2 global column numbers for first cscd
+    dn          - Number of local columns
+    dia         - First cscd starting index of each column in *ja* and *a*
+    dja         - Row of each element in first CSCD
+    da          - value of each cscd in first CSCD (can be NULL)
+    l2g         - local 2 global column numbers for first cscd
+    comm        - MPI communicator
+
+  Returns:
+    EXIT_SUCCESS if already well distributed, 2 if redistributed
+
+ */
+int z_cscd_redispatch_scotch(pastix_int_t   n, pastix_int_t *   ia, pastix_int_t *   ja, pastix_complex64_t *   a, pastix_int_t *   l2g,
+                           pastix_int_t *dn, pastix_int_t ** dia, pastix_int_t ** dja, pastix_complex64_t ** da, pastix_int_t ** dl2g,
+                           MPI_Comm comm);
+#endif /* CSCD_UTILS_INTERN_H */

z_spm.c

+/**
+ *
+ * @file z_spm.c
+ *
+ *  PaStiX spm routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @date 2013-06-24
+ *
+ * @precisions normal z -> c d s p
+ *
+ **/
+#include "common.h"
+#include "csc.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_spm
+ *
+ * z_spmSort - 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[in,out] 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( pastix_csc_t *csc )
+{
+    pastix_int_t       *colptr = csc->colptr;
+    pastix_int_t       *rowptr = csc->rowptr;
+    pastix_complex64_t *values = csc->values;
+    void *sortptr[2];
+    pastix_int_t n = csc->n;
+    pastix_int_t i, size;
+    (void)sortptr;
+
+    if (csc->dof > 1){
+        fprintf(stderr, "z_spmSort: Number of dof (%d) greater than one not supported\n", (int)csc->dof);
+        exit(1);
+    }
+
+    /* Sort in place each subset */
+    if ( csc->fmttype == PastixCSC ) {
+        for (i=0; i<n; i++, colptr++)
+        {
+            size = colptr[1] - colptr[0];
+
+#if defined(PRECISION_p)
+            intSort1asc1( rowptr, size );
+#else
+            sortptr[0] = rowptr;
+            sortptr[1] = values;
+            z_qsortIntFloatAsc( sortptr, size );
+#endif
+            rowptr += size;
+            values += size;
+        }
+    }
+    else if ( csc->fmttype == PastixCSR ) {
+        for (i=0; i<n; i++, rowptr++)
+        {
+            size = rowptr[1] - rowptr[0];
+
+#if defined(PRECISION_p)
+            intSort1asc1( colptr, size );
+#else
+            sortptr[0] = colptr;
+            sortptr[1] = values;
+            z_qsortIntFloatAsc( sortptr, size );
+#endif
+            colptr += size;
+            values += size;
+        }
+    }
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_spm
+ * @ingroup pastix_internal
+ *
+ * z_spmMergeDuplicate - 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[in,out] 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
+ *          \retval The number of vertices that were merged
+ *
+ *******************************************************************************/
+pastix_int_t
+z_spmMergeDuplicate( pastix_csc_t *csc )
+{
+    pastix_int_t       *colptr = csc->colptr;
+    pastix_int_t       *oldrow = csc->rowptr;
+    pastix_int_t       *newrow = csc->rowptr;
+    pastix_complex64_t *newval = csc->values;
+    pastix_complex64_t *oldval = csc->values;
+    pastix_int_t n       = csc->n;
+    pastix_int_t baseval = csc->colptr[0];
+    pastix_int_t dof2    = csc->dof * csc->dof;
+    pastix_int_t idx, i, j, d, size;
+    pastix_int_t merge = 0;
+    (void)d;
+
+    if ( csc->fmttype == PastixCSC ) {
+        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(pastix_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) && (csc->nnz         == idx)) ||
+                ((merge != 0) && (csc->nnz - merge == idx)) );
+
+        if (merge > 0) {
+            csc->nnz = csc->nnz - merge;
+            csc->gnnz = csc->nnz;
+
+            newrow = malloc( csc->nnz * sizeof(pastix_int_t) );
+            memcpy( newrow, csc->rowptr, csc->nnz * sizeof(pastix_int_t) );
+            free(csc->rowptr);
+            csc->rowptr = newrow;
+
+#if !defined(PRECISION_p)
+            newval = malloc( csc->nnz * dof2 * sizeof(pastix_int_t) );
+            memcpy( newval, csc->values, csc->nnz * dof2 * sizeof(pastix_complex64_t) );
+            free(csc->values);
+            csc->values = newval;
+#endif
+        }
+    }
+
+    return merge;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_spm
+ * @ingroup pastix_internal
+ *
+ * z_spmSymmetrize - This routine corrects the sparse matrix structure if it's
+ * pattern is not symmetric. It returns the new symmetric pattern with zeores on
+ * the new entries.
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          On entry, the pointer to the sparse matrix structure.
+ *          On exit, the same sparse matrix with extra entires that makes it
+ *          pattern symmetric.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *          \retval Returns the number of elements added to the matrix.
+ *
+ *******************************************************************************/
+pastix_int_t
+z_spmSymmetrize( pastix_csc_t *csc )
+{
+    pastix_complex64_t *oldval, *valtmp, *newval = NULL;
+    pastix_int_t *oldcol, *coltmp, *newcol = NULL;
+    pastix_int_t *oldrow, *rowtmp, *newrow = NULL;
+    pastix_int_t *toaddtab, *toaddtmp, toaddcnt, toaddsze;
+    pastix_int_t  n       = csc->n;
+    pastix_int_t  dof2    = csc->dof * csc->dof;
+    pastix_int_t  i, j, k, r, size;
+    pastix_int_t  baseval;
+
+    toaddcnt = 0;
+    toaddsze = 0;
+
+    if ( (csc->fmttype == PastixCSC) || (csc->fmttype == PastixCSR) ) {
+        if (csc->fmttype == PastixCSC) {
+            oldcol = csc->colptr;
+            coltmp = csc->colptr;
+            oldrow = csc->rowptr;
+            rowtmp = csc->rowptr;
+        }
+        else {
+            oldcol = csc->rowptr;
+            coltmp = csc->rowptr;
+            oldrow = csc->colptr;
+            rowtmp = csc->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) */
+                    pastix_int_t frow = oldcol[ j ] - baseval;
+                    pastix_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 csc is sorted so we will not find it later */
+                            break;
+                        }
+                    }
+
+                    if ( !found ) {
+                        if ( newcol == NULL ) {
+                            newcol = malloc( (csc->n+1) * sizeof(pastix_int_t) );
+                            for (k=0; k<csc->n; k++) {
+                                newcol[k] = oldcol[k+1] - oldcol[k];
+                            }
+
+                            /* Let's start with a buffer that will contain 5% of extra elements */
+                            toaddsze = pastix_imax( 0.05 * (double)csc->nnz, 1 );
+                            MALLOC_INTERN(toaddtab, 2*toaddsze, pastix_int_t);
+                        }
+
+                        if (toaddcnt >= toaddsze) {
+                            toaddsze *= 2;
+                            toaddtab = (pastix_int_t*)memRealloc(toaddtab, 2*toaddsze*sizeof(pastix_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) {
+            pastix_int_t newsize, oldsize;
+
+            /* Sort the array per column */
+            intSort2asc1(toaddtab, toaddcnt);
+
+            csc->nnz = csc->nnz + toaddcnt;
+            csc->gnnz = csc->nnz;
+
+            newrow = malloc( csc->nnz        * sizeof(pastix_int_t) );
+#if !defined(PRECISION_p)
+            newval = malloc( csc->nnz * dof2 * sizeof(pastix_complex64_t) );
+#endif
+            coltmp = newcol;
+            rowtmp = newrow;
+            valtmp = newval;
+            oldval = csc->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(pastix_int_t) );
+#if !defined(PRECISION_p)
+                memcpy( valtmp, oldval, oldsize * dof2 * sizeof(pastix_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(pastix_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 == csc->nnz );
+            free( toaddtab );
+            free( csc->colptr );
+            free( csc->rowptr );
+            free( csc->values );
+            if (csc->fmttype == PastixCSC) {
+                csc->colptr = newcol;
+                csc->rowptr = newrow;
+            }
+            else {
+                csc->colptr = newrow;
+                csc->rowptr = newcol;
+            }
+            csc->values = newval;
+        }
+    }
+
+    return toaddcnt;
+}
+

z_spm.h

+/**
+ *
+ * @file z_spm.h
+ *
+ *  PaStiX sparse matrix routines to handle different format of sparse matrices.
+ *  $COPYRIGHTS$
+ *
+ * @version 5.1.0
+ * @author Xavier Lacoste
+ * @author Theophile Terraz
+ * @author Pierre Ramet
+ * @author Mathieu Faverge
+ * @date 2013-06-24
+ *
+ * @precisions normal z -> c d s p
+ *
+ **/
+#ifndef _z_spm_H_
+#define _z_spm_H_
+
+int z_spmConvertCSC2CSR( pastix_csc_t *spm );
+int z_spmConvertCSC2IJV( pastix_csc_t *spm );
+int z_spmConvertCSR2CSC( pastix_csc_t *spm );
+int z_spmConvertCSR2IJV( pastix_csc_t *spm );
+int z_spmConvertIJV2CSC( pastix_csc_t *spm );
+int z_spmConvertIJV2CSR( pastix_csc_t *spm );
+
+int    z_spmGeCSCv(int trans, pastix_complex64_t alpha, const pastix_csc_t *csc, const pastix_complex64_t *x, pastix_complex64_t beta, pastix_complex64_t *y);
+int    z_spmSyCSCv(           pastix_complex64_t alpha, const pastix_csc_t *csc, const pastix_complex64_t *x, pastix_complex64_t beta, pastix_complex64_t *y);
+int    z_spmHeCSCv(           pastix_complex64_t alpha, const pastix_csc_t *csc, const pastix_complex64_t *x, pastix_complex64_t beta, pastix_complex64_t *y);
+
+double       z_spmNorm( int ntype, const pastix_csc_t *csc );
+void         z_spmSort( pastix_csc_t *csc );
+pastix_int_t z_spmMergeDuplicate( pastix_csc_t *csc );
+pastix_int_t z_spmSymmetrize( pastix_csc_t *csc );
+
+int z_spmGenRHS(int type, int nrhs, const pastix_csc_t *spm, void *x, int ldx, void *b, int ldb );
+int z_spmCheckAxb( int nrhs, const pastix_csc_t *spm, void *x0, int ldx0, void *b, int ldb, const void *x, int ldx );
+
+#endif /* _z_spm_H_ */

z_spm_convert_to_csc.c

+/**
+ *
+ * @file z_spm_convert_to_csc.c
+ *
+ *  PaStiX csc routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @date 2015-01-01
+ *
+ * @precisions normal z -> c d s p
+ **/
+#include "common.h"
+#include "csc.h"
+#include "z_spm.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc
+ *
+ * z_spmConvertIJV2CSC - convert a matrix in IJV format to a matrix in CSC
+ * format.
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          The ijv matrix at enter,
+ *          the csc matrix at exit.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS
+ *
+ *******************************************************************************/
+int
+z_spmConvertIJV2CSC( pastix_csc_t *spm )
+{
+#if !defined(PRECISION_p)
+    pastix_complex64_t *navals = NULL;
+    pastix_complex64_t *oavals = NULL;
+#endif
+    pastix_int_t       *spmptx, *otmp;
+    pastix_int_t i, j, tmp, baseval, total;
+    pastix_csc_t oldspm;
+
+    /* Backup the input */
+    memcpy( &oldspm, spm, sizeof(pastix_csc_t) );
+
+    /*
+     * Check the baseval, we consider that arrays are sorted by columns or rows
+     */
+    baseval = spmFindBase( spm );
+
+    /* Compute the new colptr */
+    spm->colptr = (pastix_int_t *) calloc(spm->n+1,sizeof(pastix_int_t));
+
+    /* Compute the number of edges per row */
+    spmptx = spm->colptr - baseval;
+    otmp   = oldspm.colptr;
+    for (i=0; i<spm->nnz; i++, otmp++)
+    {
+        spmptx[ *otmp ] ++;
+    }
+
+    /* Compute the indexes in C numbering for the following sort */
+    total = 0;
+    spmptx = spm->colptr;
+    for (i=0; i<(spm->n+1); i++, spmptx++)
+    {
+        tmp = *spmptx;
+        *spmptx = total;
+        total += tmp;
+    }
+    assert( total == spm->nnz );
+
+    /* Sort the rows and avals arrays by column */
+    spm->rowptr = malloc(spm->nnz * sizeof(pastix_int_t));
+
+#if defined(PRECISION_p)
+    spm->values = NULL;
+#else
+    spm->values = malloc(spm->nnz * sizeof(pastix_complex64_t));
+    navals = (pastix_complex64_t*)(spm->values);
+    oavals = (pastix_complex64_t*)(oldspm.values);
+#endif
+
+    for (j=0; j<spm->nnz; j++)
+    {
+        i = oldspm.colptr[j] - baseval;
+
+        spm->rowptr[ spm->colptr[i] ] = oldspm.rowptr[j];
+
+#if !defined(PRECISION_p)
+        navals[ spm->colptr[i] ] = oavals[j];
+#endif
+        (spm->colptr[i])++;
+
+        assert( spm->colptr[i] <= spm->colptr[i+1] );
+    }
+
+    /* Rebuild the colptr with the correct baseval */
+    tmp = spm->colptr[0];
+    spm->colptr[0] = baseval;
+
+    spmptx = spm->colptr + 1;
+    for (i=1; i<(spm->n+1); i++, spmptx++)
+    {
+        total = *spmptx;
+        *spmptx = tmp + baseval;
+        tmp = total;
+    }
+    assert( spm->colptr[ spm->n ] == (spm->nnz+baseval) );
+
+    spmExit( &oldspm );
+
+    spm->fmttype = PastixCSC;
+
+    return PASTIX_SUCCESS;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc
+ *
+ * z_spmConvertCSR2CSC - convert a matrix in CSR format to a matrix in CSC
+ * format. If the matrix is PastixSymmetric or PastixHermitian, then the
+ * transpose or respectively the conjugate is returned.
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          The csr matrix at enter,
+ *          the csc matrix at exit.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS
+ *
+ *******************************************************************************/
+int
+z_spmConvertCSR2CSC( pastix_csc_t *spm )
+{
+    assert( spm->loc2glob == NULL );
+    assert( spm->fmttype == PastixCSR );
+
+    spm->fmttype = PastixCSC;
+
+    switch( spm->mtxtype ) {
+#if defined(PRECISION_z) || defined(PRECISION_c)
+    case PastixHermitian:
+    {
+        /* Similar to PastixSymmetric case with conjugate of the values */
+        pastix_complex64_t *valptr = spm->values;
+        pastix_int_t i;
+
+        for(i=0; i<spm->nnz; i++, valptr++){
+            *valptr = conj( *valptr );
+        }
+    }
+#endif
+    case PastixSymmetric:
+    {
+        pastix_int_t *tmp;
+
+        /* Just need to swap the pointers */
+        tmp          = spm->rowptr;
+        spm->rowptr  = spm->colptr;
+        spm->colptr  = tmp;
+        spm->fmttype = PastixCSC;
+
+        return PASTIX_SUCCESS;
+    }
+    break;
+
+    case PastixGeneral:
+    default:
+    {
+        pastix_int_t       *row_csc;
+        pastix_int_t       *col_csc;
+#if !defined(PRECISION_p)
+        pastix_complex64_t *val_csc;
+        pastix_complex64_t *valptr = (pastix_complex64_t*)(spm->values);
+#endif
+        pastix_int_t j, k, col, row, nnz, baseval;
+
+        baseval = spmFindBase( spm );
+        nnz = spm->nnz;
+
+        row_csc = malloc(nnz * sizeof(pastix_int_t));
+        col_csc = calloc(spm->n+1,sizeof(pastix_int_t));
+
+        assert( row_csc );
+        assert( col_csc );
+
+#if !defined(PRECISION_p)
+        val_csc = malloc(nnz*sizeof(pastix_complex64_t));
+        assert( val_csc );
+#endif
+
+        /* Count the number of elements per column */
+        for (j=0; j<nnz; j++) {
+            col = spm->colptr[j] - baseval;
+            assert(col < spm->n );
+            col_csc[ col+1 ] ++;
+        }
+
+        /* Compute the index of each column */
+        col_csc[0] = 0;
+        for (j=0; j<spm->n; j++){
+            col_csc[j+1] += col_csc[j];
+        }
+
+        assert( (col_csc[spm->gN]) == nnz );
+
+        for (row=0; row<spm->n; row++) {
+            pastix_int_t fcol = spm->rowptr[row  ] - baseval;
+            pastix_int_t lcol = spm->rowptr[row+1] - baseval;
+
+            for (k=fcol; k<lcol; k++) {
+                col = spm->colptr[k] - baseval;
+                j = col_csc[col];
+                row_csc[j] = row + baseval;
+
+#if !defined(PRECISION_p)
+                val_csc[j] = valptr[k];
+#endif
+                col_csc[col] ++;
+            }
+        }
+
+        /* Restore the colptr indexes */
+        {
+            pastix_int_t tmp, tmp2;
+
+            tmp = col_csc[0];
+            col_csc[0] = baseval;
+            for (j=0; j<spm->n; j++) {
+                tmp2 = col_csc[j+1];
+                col_csc[j+1] = tmp + baseval;
+                tmp = tmp2;
+            }
+        }
+
+        spmExit( spm );
+        spm->colptr = col_csc;
+        spm->rowptr = row_csc;
+#if !defined(PRECISION_p)
+        spm->values = val_csc;
+#else
+        spm->values = NULL;
+#endif
+    }
+    }
+
+    return PASTIX_SUCCESS;
+}

z_spm_convert_to_csr.c

+/**
+ *
+ * @file z_spm_convert_to_csr.c
+ *
+ *  PaStiX csc routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @date 2015-01-01
+ *
+ * @precisions normal z -> c d s p
+ **/
+#include "common.h"
+#include "csc.h"
+#include "z_spm.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc
+ *
+ * z_spmConvertCSC2CSR - convert a matrix in CSC format to a matrix in CSR
+ * format. If the matrix is PastixSymmetric or PastixHermitian, then the
+ * transpose or respectively the conjugate is returned.
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          The csc matrix at enter,
+ *          the csr matrix at exit.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS
+ *
+ *******************************************************************************/
+int
+z_spmConvertCSC2CSR( pastix_csc_t *spm )
+{
+    pastix_int_t *tmp;
+    pastix_int_t  result;
+
+    switch( spm->mtxtype ) {
+#if defined(PRECISION_z) || defined(PRECISION_c)
+    case PastixHermitian:
+    {
+        /* Similar to PastixSymmetric case with conjugate of the values */
+        pastix_complex64_t *valptr = spm->values;
+        pastix_int_t i;
+
+        for(i=0; i<spm->nnz; i++, valptr++){
+            *valptr = conj( *valptr );
+        }
+    }
+#endif
+    case PastixSymmetric:
+    {
+        pastix_int_t *tmp;
+
+        /* Just need to swap the pointers */
+        tmp          = spm->rowptr;
+        spm->rowptr  = spm->colptr;
+        spm->colptr  = tmp;
+        spm->fmttype = PastixCSR;
+
+        return PASTIX_SUCCESS;
+    }
+    break;
+
+    case PastixGeneral:
+    default:
+    {
+        /* Transpose the spm in CSC to trans(spm) in CSR */
+        tmp          = spm->rowptr;
+        spm->rowptr  = spm->colptr;
+        spm->colptr  = tmp;
+        spm->fmttype = PastixCSR;
+
+        /* Convert trans(spm) in CSR to trans(spm) in CSC */
+        result = z_spmConvertCSR2CSC( spm );
+
+        /* Transpose trans(spm) in CSC to obtain the spm in CSR */
+        tmp          = spm->rowptr;
+        spm->rowptr  = spm->colptr;
+        spm->colptr  = tmp;
+        spm->fmttype = PastixCSR;
+    }
+    }
+
+    return result;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc
+ *
+ * z_spmConvertIJV2CSR - convert a matrix in IJV format to a matrix in CSR
+ * format.
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          The ijv matrix at enter,
+ *          the csr matrix at exit.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS
+ *
+ *******************************************************************************/
+int
+z_spmConvertIJV2CSR( pastix_csc_t *spm )
+{
+#if !defined(PRECISION_p)
+    pastix_complex64_t *navals = NULL;
+    pastix_complex64_t *oavals = NULL;
+#endif
+    pastix_int_t       *spmptx, *otmp;
+    pastix_int_t i, j, tmp, baseval, total;
+    pastix_csc_t oldspm;
+
+    /* Backup the input */
+    memcpy( &oldspm, spm, sizeof(pastix_csc_t) );
+
+    /*
+     * Check the baseval, we consider that arrays are sorted by columns or rows
+     */
+    baseval = spmFindBase( spm );
+
+    /* Compute the new rowptr */
+    spm->rowptr = (pastix_int_t *) calloc(spm->n+1,sizeof(pastix_int_t));
+
+    /* Compute the number of edges per row */
+    spmptx = spm->rowptr - baseval;
+    otmp   = oldspm.rowptr;
+    for (i=0; i<spm->nnz; i++, otmp++)
+    {
+        spmptx[ *otmp ] ++;
+    }
+
+    /* Compute the indexes in C numbering for the following sort */
+    total = 0;
+    spmptx = spm->rowptr;
+    for (i=0; i<(spm->n+1); i++, spmptx++)
+    {
+        tmp = *spmptx;
+        *spmptx = total;
+        total += tmp;
+    }
+    assert( total == spm->nnz );
+
+    /* Sort the colptr and avals arrays by rows */
+    spm->colptr  = malloc(spm->nnz * sizeof(pastix_int_t));
+
+#if defined(PRECISION_p)
+    spm->values = NULL;
+#else
+    spm->values = malloc(spm->nnz * sizeof(pastix_complex64_t));
+    navals = (pastix_complex64_t*)(spm->values);
+    oavals = (pastix_complex64_t*)(oldspm.values);
+#endif
+
+    for (j=0; j<spm->nnz; j++)
+    {
+        i = oldspm.rowptr[j] - baseval;
+
+        spm->colptr[ spm->rowptr[i] ] = oldspm.colptr[j];
+
+#if !defined(PRECISION_p)
+        navals[ spm->rowptr[i] ] = oavals[j];
+#endif
+        (spm->rowptr[i])++;
+
+        assert( spm->rowptr[i] <= spm->rowptr[i+1] );
+    }
+
+    /* Rebuild the rows (rowptr) with the correct baseval */
+    tmp = spm->rowptr[0];
+    spm->rowptr[0] = baseval;
+
+    spmptx = spm->rowptr + 1;
+    for (i=1; i<(spm->n+1); i++, spmptx++)
+    {
+        total = *spmptx;
+        *spmptx = tmp + baseval;
+        tmp = total;
+    }
+    assert( spm->rowptr[ spm->n ] == (spm->nnz+baseval) );
+
+    spmExit( &oldspm );
+
+    spm->fmttype = PastixCSR;
+
+    return PASTIX_SUCCESS;
+}

z_spm_convert_to_ijv.c

+/**
+ *
+ * @file z_spm_convert_to_ijv.c
+ *
+ *  PaStiX csc routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @date 2015-01-01
+ *
+ * @precisions normal z -> c d s p
+ **/
+#include "common.h"
+#include "csc.h"
+#include "z_spm.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc
+ *
+ * z_spmConvertCSC2IJV - convert a matrix in CSC format to a matrix in IJV
+ * format.
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          The csc matrix at enter,
+ *          the ijv matrix at exit.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS
+ *
+ *******************************************************************************/
+int
+z_spmConvertCSC2IJV( pastix_csc_t *spm )
+{
+    pastix_int_t *col_ijv, *colptr;
+    pastix_int_t i, j, baseval, nnz;
+
+    /*
+     * Check the baseval
+     */
+    baseval = spmFindBase( spm );
+    nnz = spm->nnz;
+    spm->fmttype = PastixIJV;
+
+    col_ijv = malloc(nnz*sizeof(pastix_int_t));
+    assert( col_ijv );
+
+    colptr = col_ijv;
+    for(i=0; i<spm->n; i++)
+    {
+        for(j=spm->colptr[i]; j<spm->colptr[i+1]; j++)
+        {
+            *colptr = i+baseval; colptr++;
+        }
+    }
+
+    memFree_null(spm->colptr);
+    spm->colptr = col_ijv;
+
+    return PASTIX_SUCCESS;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc
+ *
+ * z_spmConvertCSR2IJV - convert a matrix in CSR format to a matrix in IJV
+ * format.
+ *
+ *******************************************************************************
+ *
+ * @param[in,out] spm
+ *          The csr matrix at enter,
+ *          the ijv matrix at exit.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS
+ *
+ *******************************************************************************/
+int
+z_spmConvertCSR2IJV( pastix_csc_t *spm )
+{
+    pastix_int_t *row_ijv, *rowptr;
+    pastix_int_t i, j, baseval, nnz;
+
+    /*
+     * Check the baseval
+     */
+    baseval = spmFindBase( spm );
+    nnz = spm->nnz;
+    spm->fmttype = PastixIJV;
+
+    row_ijv = malloc(nnz*sizeof(pastix_int_t));
+    assert( row_ijv );
+
+    rowptr = row_ijv;
+    for(i=0; i<spm->n; i++)
+    {
+        for(j=spm->rowptr[i]; j<spm->rowptr[i+1]; j++)
+        {
+            *rowptr = i+baseval; rowptr++;
+        }
+    }
+
+    memFree_null(spm->rowptr);
+    spm->rowptr = row_ijv;
+
+    return PASTIX_SUCCESS;
+}