Merge default and integrate fixes

CMakeLists.txt

@@ -4,6 +4,8 @@ include(RulesPrecisions)
 set(generated_sources "")
 set(generated_headers "")
 
+include_directories("${CMAKE_CURRENT_SOURCE_DIR}/drivers")
+
 set(HEADERS
   z_spm.h
 )
@@ -15,17 +17,25 @@ precisions_rules_py(generated_headers
 
 set(SOURCES
   z_spm.c
+  z_spm_2dense.c
+  z_spm_dof_extend.c
+  z_spm_norm.c
+
   z_spm_convert_to_csc.c
   z_spm_convert_to_csr.c
   z_spm_convert_to_ijv.c
+  z_spm_expand.c
   z_spm_genrhs.c
+  z_spm_integer.c
+  z_spm_laplacian.c
   z_spm_matrixvector.c
-  z_spm_norm.c
+  z_spm_print.c
   )
 
 set(spm_headers
   ${generated_headers}
   spm.h
+  spm_drivers.h
   )
 
 ### Generate the sources in all precisions
@@ -37,17 +47,35 @@ set(spm_sources
   ${generated_sources}
   spm.c
   spm_io.c
+  spm_integers.c
+  spm_dofs.c
+  spm_read_driver.c
+  drivers/skitf.f
+  drivers/iohb.c
+  drivers/mmio.c
+  drivers/laplacian.c
+  drivers/readhb.c
+  drivers/readijv.c
+  drivers/readmm.c
+  drivers/readrsa.c
   )
 
 add_custom_target(spm_headers_tgt
   DEPENDS ${spm_headers} )
 
-set_source_files_properties(
-  spm.c
-  PROPERTIES DEPENDS spm_headers_tgt
-  )
-
-set(PASTIX_LIB_SRCS
-  ${PASTIX_LIB_SRCS}
+add_library(pastix_spm
   ${spm_sources}
   )
+
+add_dependencies(pastix_spm
+  spm_headers_tgt
+)
+
+### Generate the lib
+if (MPI_FOUND)
+  set_target_properties(pastix_spm PROPERTIES COMPILE_FLAGS "${MPI_COMPILE_FLAGS}")
+endif (MPI_FOUND)
+
+install(TARGETS pastix_spm
+  ARCHIVE DESTINATION lib
+  LIBRARY DESTINATION lib)

drivers/iohb.h

+#ifndef IOHB_H
+#define IOHB_H
+
+int readHB_info(const char* filename, int* M, int* N, int* nz, char** Type,
+                                                      int* Nrhs);
+
+int readHB_header(FILE* in_file, char* Title, char* Key, char* Type,
+                    int* Nrow, int* Ncol, int* Nnzero, int* Nrhs,
+                    char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
+                    int* Ptrcrd, int* Indcrd, int* Valcrd, int* Rhscrd,
+                    char *Rhstype);
+
+int readHB_mat_double(const char* filename, int colptr[], int rowind[],
+                                                                 double val[]);
+
+int readHB_newmat_double(const char* filename, int* M, int* N, int* nonzeros,
+                         int** colptr, int** rowind, double** val);
+
+int readHB_aux_double(const char* filename, const char AuxType, double b[]);
+
+int readHB_newaux_double(const char* filename, const char AuxType, double** b);
+
+int writeHB_mat_double(const char* filename, int M, int N,
+                        int nz, const int colptr[], const int rowind[],
+                        const double val[], int Nrhs, const double rhs[],
+                        const double guess[], const double exact[],
+                        const char* Title, const char* Key, const char* Type,
+                        char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
+                        const char* Rhstype);
+
+int readHB_mat_char(const char* filename, int colptr[], int rowind[],
+                                           char val[], char* Valfmt);
+
+int readHB_newmat_char(const char* filename, int* M, int* N, int* nonzeros, int** colptr,
+                          int** rowind, char** val, char** Valfmt);
+
+int readHB_aux_char(const char* filename, const char AuxType, char b[]);
+
+int readHB_newaux_char(const char* filename, const char AuxType, char** b, char** Rhsfmt);
+
+int writeHB_mat_char(const char* filename, int M, int N,
+                        int nz, const int colptr[], const int rowind[],
+                        const char val[], int Nrhs, const char rhs[],
+                        const char guess[], const char exact[],
+                        const char* Title, const char* Key, const char* Type,
+                        char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
+                        const char* Rhstype);
+
+int ParseIfmt(char* fmt, int* perline, int* width);
+
+int ParseRfmt(char* fmt, int* perline, int* width, int* prec, char* flag);
+
+void IOHBTerminate(char* message);
+
+#endif

drivers/laplacian.c

+/**
+ * @file laplacian.c
+ *
+ *  This file contains the user routine to generate 1, 2 and 3 dimensional Laplacian
+ *  matrices.
+ *  $COPYRIGHTS$
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Pierre Ramet
+ * @author Xavier Lacoste
+ * @author Theophile Terraz
+ * @date 2011-11-11
+ *
+ **/
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "common.h"
+#include "spm_drivers.h"
+#include "drivers/laplacian.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * laplacian_usage - Print the usage information to generate correct Laplacian
+ * matrices.
+ *
+ *******************************************************************************/
+static inline void
+laplacian_usage(void)
+{
+    fprintf(stderr,
+            "Usage: genLaplacian( \"[<type>:]<dim1>[:<dim2>[:<dim3>]]\" )\n"
+            "   <type> p = pattern only\n"
+            "          s = real simple\n"
+            "          d = real double [default]\n"
+            "          c = complex simple\n"
+            "          z = complex double\n"
+            "   <dim1> size of the first dimension of the 1D|2D|3D laplacian\n"
+            "   <dim2> size of the second dimension of the 2D|3D laplacian\n"
+            "   <dim3> size of the third dimension of the 3D laplacian\n"
+            "   Example:\n"
+            "     genLaplacian( \"z:10:20\" ) generates a 2D complex laplacian matrix of size 200.\n"
+            );
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * laplacian_parse_info - Parse information given through the filename string to
+ * configure the laplacian matrix to generate.
+ *
+ *******************************************************************************
+ *
+ * @param[in] filename
+ *          Configuration string of the Laplacian. See laplacian_usage() for
+ *          more information.
+ *
+ * @param[in,out] csc
+ *          At start, an allocated csc structure that will store the Lapalcian
+ *          matrix.
+ *          At exit, the fields of the csc are initialized and especially the
+ *          type, symmetry and number of unknows are setup.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been generated successfully
+ *      \retval PASTIX_ERR_BADPARAMETER if the configuration string is incorrect
+ *
+ *******************************************************************************/
+static inline int
+laplacian_parse_info( const char   *filename,
+                      pastix_csc_t *csc,
+                      pastix_int_t *dim1,
+                      pastix_int_t *dim2,
+                      pastix_int_t *dim3 )
+{
+    long tmp1, tmp2, tmp3;
+    csc->colptr   = NULL;
+    csc->rowptr   = NULL;
+    csc->values   = NULL;
+    csc->loc2glob = NULL;
+
+    /* Look for the datatype */
+    {
+        char flt;
+        char *tmpf = strdup( filename );
+
+        if ( sscanf( filename, "%c:%s", &flt, tmpf ) == 2 ) {
+            filename += 2;
+            switch( flt ){
+            case 'Z':
+            case 'z':
+                csc->flttype = PastixComplex64;
+                break;
+
+            case 'C':
+            case 'c':
+                csc->flttype = PastixComplex32;
+                break;
+
+            case 'D':
+            case 'd':
+                csc->flttype = PastixDouble;
+                break;
+
+            case 'S':
+            case 's':
+                csc->flttype = PastixFloat;
+                break;
+
+            case 'P':
+            case 'p':
+                csc->flttype = PastixPattern;
+                break;
+
+            case '1':
+            case '2':
+            case '3':
+            case '4':
+            case '5':
+            case '6':
+            case '7':
+            case '8':
+            case '9':
+                csc->flttype = PastixDouble;
+                /*
+                 * The first dimension is only one character long so we come
+                 * back to the beginning of the string
+                 */
+                filename -= 2;
+                break;
+
+            default:
+                laplacian_usage();
+                return PASTIX_ERR_BADPARAMETER;
+            }
+        }
+        else {
+            csc->flttype = PastixDouble;
+        }
+
+        free(tmpf);
+    }
+
+    /* Scan the dimensions */
+    *dim1 = *dim2 = *dim3 = 0;
+
+    if ( sscanf( filename, "%ld:%ld:%ld", &tmp1, &tmp2, &tmp3 ) == 3 ) {
+        *dim1 = (pastix_int_t)tmp1;
+        *dim2 = (pastix_int_t)tmp2;
+        *dim3 = (pastix_int_t)tmp3;
+        csc->gN = (*dim1)*(*dim2)*(*dim3);
+    }
+    else if ( sscanf( filename, "%ld:%ld", &tmp1, &tmp2 ) == 2 ) {
+        *dim1 = (pastix_int_t)tmp1;
+        *dim2 = (pastix_int_t)tmp2;
+        csc->gN = (*dim1)*(*dim2);
+    }
+    else if ( sscanf( filename, "%ld", &tmp1 ) == 1 ) {
+        *dim1 = (pastix_int_t)tmp1;
+        csc->gN = *dim1;
+    }
+    else {
+        laplacian_usage();
+        return PASTIX_ERR_BADPARAMETER;
+    }
+
+    /* One of the dimension was set to 0 */
+    if ( csc->gN == 0 ) {
+        laplacian_usage();
+        return PASTIX_ERR_BADPARAMETER;
+    }
+
+    csc->n = csc->gN;
+    return PASTIX_SUCCESS;
+}
+
+static void (*laplacian_table1D[6])(pastix_csc_t*, pastix_int_t) =
+{
+    p_spmLaplacian1D,
+    NULL,
+    s_spmLaplacian1D,
+    d_spmLaplacian1D,
+    c_spmLaplacian1D,
+    z_spmLaplacian1D
+};
+
+static void (*laplacian_table2D[6])(pastix_csc_t*, pastix_int_t, pastix_int_t) =
+{
+    p_spmLaplacian2D,
+    NULL,
+    s_spmLaplacian2D,
+    d_spmLaplacian2D,
+    c_spmLaplacian2D,
+    z_spmLaplacian2D
+};
+
+static void (*laplacian_table3D[6])(pastix_csc_t*, pastix_int_t, pastix_int_t, pastix_int_t) =
+{
+    p_spmLaplacian3D,
+    NULL,
+    s_spmLaplacian3D,
+    d_spmLaplacian3D,
+    c_spmLaplacian3D,
+    z_spmLaplacian3D
+};
+
+static void (*extended_laplacian_table2D[6])(pastix_csc_t*, pastix_int_t, pastix_int_t) =
+{
+    p_spmLaplacian2D,
+    NULL,
+    s_spmExtendedLaplacian2D,
+    d_spmExtendedLaplacian2D,
+    c_spmExtendedLaplacian2D,
+    z_spmExtendedLaplacian2D
+};
+
+static void (*extended_laplacian_table3D[6])(pastix_csc_t*, pastix_int_t, pastix_int_t, pastix_int_t) =
+{
+    p_spmExtendedLaplacian3D,
+    NULL,
+    s_spmExtendedLaplacian3D,
+    d_spmExtendedLaplacian3D,
+    c_spmExtendedLaplacian3D,
+    z_spmExtendedLaplacian3D
+};
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * genLaplacian - Generate a Laplacian of size csc->n
+ *
+ *******************************************************************************
+ *
+ * @param[in] filename
+ *          Configuration string of the Laplacian.
+ *          [<type>:]<dim1>[:<dim2>[:<dim3>]]
+ *             <type> p = pattern only\n"
+ *                    s = real simple\n"
+ *                    d = real double [default]\n"
+ *                    c = complex simple\n"
+ *                    z = complex double\n"
+ *             <dim1> size of the first dimension of the 1D|2D|3D laplacian\n"
+ *             <dim2> size of the second dimension of the 2D|3D laplacian\n"
+ *             <dim3> size of the third dimension of the 3D laplacian\n"
+ *
+ * @param[in,out] csc
+ *          At start, an allocated csc structure.
+ *          At exit, contains a laplacian matrix in the csc format.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been generated successfully
+ *      \retval PASTIX_ERR_BADPARAMETER if the configuration string is incorrect
+ *
+ *******************************************************************************/
+int
+genLaplacian( const char    *filename,
+              pastix_csc_t  *csc )
+{
+    pastix_int_t dim1, dim2, dim3;
+    int rc;
+
+    rc = laplacian_parse_info(filename, csc, &dim1, &dim2, &dim3);
+    if (rc != PASTIX_SUCCESS)
+        return rc;
+
+    if( dim3 > 0 ) {
+        laplacian_table3D[csc->flttype](csc, dim1, dim2, dim3);
+    }
+    else if (dim2 > 0) {
+        laplacian_table2D[csc->flttype](csc, dim1, dim2);
+    }
+    else {
+        laplacian_table1D[csc->flttype](csc, dim1);
+    }
+
+    return PASTIX_SUCCESS;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * genExtendedLaplacian - Generate a extended Laplacian of size csc->n
+ *
+ *******************************************************************************
+ *
+ * @param[in] filename
+ *          Configuration string of the Laplacian.
+ *          [<type>:]<dim1>[:<dim2>[:<dim3>]]
+ *             <type> p = pattern only
+ *                    s = real simple
+ *                    d = real double [default]
+ *                    c = complex simple
+ *                    z = complex double
+ *             <dim1> size of the first dimension of the 1D|2D|3D laplacian
+ *             <dim2> size of the second dimension of the 2D|3D laplacian
+ *             <dim3> size of the third dimension of the 3D laplacian
+ *
+ * @param[in,out] csc
+ *          At start, an allocated csc structure.
+ *          At exit, contains a laplacian matrix in the csc format.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been generated successfully
+ *      \retval PASTIX_ERR_BADPARAMETER if the configuration string is incorrect
+ *
+ *******************************************************************************/
+int
+genExtendedLaplacian( const char    *filename,
+                      pastix_csc_t  *csc )
+{
+    pastix_int_t dim1, dim2, dim3;
+    int rc;
+
+    rc = laplacian_parse_info(filename, csc, &dim1, &dim2, &dim3);
+    if (rc != PASTIX_SUCCESS)
+        return rc;
+
+    if( dim3 > 0 ) {
+        extended_laplacian_table3D[csc->flttype](csc, dim1, dim2, dim3);
+    }
+    else if (dim2 > 0) {
+        extended_laplacian_table2D[csc->flttype](csc, dim1, dim2);
+    }
+    else {
+        laplacian_table1D[csc->flttype](csc, dim1);
+    }
+
+    return PASTIX_SUCCESS;
+}

drivers/laplacian.h

+/**
+ * @file laplacian.h
+ *
+ *  $COPYRIGHTS$
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @date 2011-11-11
+ *
+ **/
+#ifndef _LAPLACIAN_H_
+#define _LAPLACIAN_H_
+
+void z_spmLaplacian1D( pastix_csc_t *csc, pastix_int_t dim1 );
+void c_spmLaplacian1D( pastix_csc_t *csc, pastix_int_t dim1 );
+void d_spmLaplacian1D( pastix_csc_t *csc, pastix_int_t dim1 );
+void s_spmLaplacian1D( pastix_csc_t *csc, pastix_int_t dim1 );
+void p_spmLaplacian1D( pastix_csc_t *csc, pastix_int_t dim1 );
+
+void z_spmLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+void c_spmLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+void d_spmLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+void s_spmLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+void p_spmLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+
+void z_spmLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+void c_spmLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+void d_spmLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+void s_spmLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+void p_spmLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+
+void z_spmExtendedLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+void c_spmExtendedLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+void d_spmExtendedLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+void s_spmExtendedLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+void p_spmExtendedLaplacian2D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2 );
+
+void z_spmExtendedLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+void c_spmExtendedLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+void d_spmExtendedLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+void s_spmExtendedLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+void p_spmExtendedLaplacian3D( pastix_csc_t *csc, pastix_int_t dim1, pastix_int_t dim2, pastix_int_t dim3 );
+
+#endif /* _LAPLACIAN_H_ */

drivers/mmio.c

+/*
+*   Matrix Market I/O library for ANSI C
+*
+*   See http://math.nist.gov/MatrixMarket for details.
+*
+*
+*/
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <ctype.h>
+
+#include "common.h"
+#include "spm_drivers.h"
+#include "drivers/mmio.h"
+
+int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
+                               double **val_, int **row_, int **col_)
+{
+    FILE *f;
+    MM_typecode matcode;
+    int M, N, nz;
+    int i;
+    double *val;
+    int *row, *col;
+
+    if ((f = fopen(fname, "r")) == NULL)
+            return -1;
+
+
+    if (mm_read_banner(f, &matcode) != 0)
+    {
+        printf("mm_read_unsymetric: Could not process Matrix Market banner ");
+        printf(" in file [%s]\n", fname);
+        return -1;
+    }
+
+
+
+    if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
+            mm_is_sparse(matcode)))
+    {
+        fprintf(stderr, "Sorry, this application does not support ");
+        fprintf(stderr, "Market Market type: [%s]\n",
+                mm_typecode_to_str(matcode));
+        return -1;
+    }
+
+    /* find out size of sparse matrix: M, N, nz .... */
+
+    if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
+    {
+        fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
+        return -1;
+    }
+
+    *M_ = M;
+    *N_ = N;
+    *nz_ = nz;
+
+    /* reseve memory for matrices */
+
+    row = (int *) malloc(nz * sizeof(int));
+    col = (int *) malloc(nz * sizeof(int));
+    val = (double *) malloc(nz * sizeof(double));
+
+    *val_ = val;
+    *row_ = row;
+    *col_ = col;
+
+    /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
+    /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
+    /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
+
+    for (i=0; i<nz; i++)
+    {
+      if (3 != fscanf(f, "%d %d %lg\n", &row[i], &col[i], &val[i]))
+        {
+          fprintf(stderr, "Error: reading matrix\n");
+          return 1;
+        }
+        row[i]--;  /* adjust from 1-based to 0-based */
+        col[i]--;
+    }
+    fclose(f);
+
+    return 0;
+}
+
+int mm_is_valid(MM_typecode matcode)
+{
+    if (!mm_is_matrix(matcode)) return 0;
+    if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
+    if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
+    if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) ||
+                mm_is_skew(matcode))) return 0;
+    return 1;
+}
+
+int mm_read_banner(FILE *f, MM_typecode *matcode)
+{
+    char line[MM_MAX_LINE_LENGTH];
+    char banner[MM_MAX_TOKEN_LENGTH];
+    char mtx[MM_MAX_TOKEN_LENGTH];
+    char crd[MM_MAX_TOKEN_LENGTH];
+    char data_type[MM_MAX_TOKEN_LENGTH];
+    char storage_scheme[MM_MAX_TOKEN_LENGTH];
+    char *p;
+
+
+    mm_clear_typecode(matcode);
+
+    if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
+        return MM_PREMATURE_EOF;
+
+    if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type,
+        storage_scheme) != 5)
+        return MM_PREMATURE_EOF;
+
+    for (p=mtx; *p!='\0'; *p=(char)tolower((int)*p),p++);  /* convert to lower case */
+    for (p=crd; *p!='\0'; *p=(char)tolower((int)*p),p++);
+    for (p=data_type; *p!='\0'; *p=(char)tolower((int)*p),p++);
+    for (p=storage_scheme; *p!='\0'; *p=(char)tolower((int)*p),p++);
+
+    /* check for banner */
+    if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
+        return MM_NO_HEADER;
+
+    /* first field should be "mtx" */
+    if (strcmp(mtx, MM_MTX_STR) != 0)
+        return  MM_UNSUPPORTED_TYPE;
+    mm_set_matrix(matcode);
+
+
+    /* second field describes whether this is a sparse matrix (in coordinate
+            storgae) or a dense array */
+
+
+    if (strcmp(crd, MM_SPARSE_STR) == 0)
+        mm_set_sparse(matcode);
+    else
+    if (strcmp(crd, MM_DENSE_STR) == 0)
+            mm_set_dense(matcode);
+    else
+        return MM_UNSUPPORTED_TYPE;
+
+
+    /* third field */
+
+    if (strcmp(data_type, MM_REAL_STR) == 0)
+        mm_set_real(matcode);
+    else
+    if (strcmp(data_type, MM_COMPLEX_STR) == 0)
+        mm_set_complex(matcode);
+    else
+    if (strcmp(data_type, MM_PATTERN_STR) == 0)
+        mm_set_pattern(matcode);
+    else
+    if (strcmp(data_type, MM_INT_STR) == 0)
+        mm_set_integer(matcode);
+    else
+        return MM_UNSUPPORTED_TYPE;
+
+
+    /* fourth field */
+
+    if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
+        mm_set_general(matcode);
+    else
+    if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
+        mm_set_symmetric(matcode);
+    else
+    if (strcmp(storage_scheme, MM_HERM_STR) == 0)
+        mm_set_hermitian(matcode);
+    else
+    if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
+        mm_set_skew(matcode);
+    else
+        return MM_UNSUPPORTED_TYPE;
+
+
+    return 0;
+}
+
+int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
+{
+    if (fprintf(f, "%d %d %d\n", M, N, nz) != 3)
+        return MM_COULD_NOT_WRITE_FILE;
+    else
+        return 0;
+}
+
+int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
+{
+    char line[MM_MAX_LINE_LENGTH];
+    int num_items_read;
+
+    /* set return null parameter values, in case we exit with errors */
+    *M = *N = *nz = 0;
+
+    /* now continue scanning until you reach the end-of-comments */
+    do
+    {
+        if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
+            return MM_PREMATURE_EOF;
+    }while (line[0] == '%');
+
+    /* line[] is either blank or has M,N, nz */
+    if (sscanf(line, "%d %d %d", M, N, nz) == 3)
+        return 0;
+
+    else
+    do
+    {
+        num_items_read = fscanf(f, "%d %d %d", M, N, nz);
+        if (num_items_read == EOF) return MM_PREMATURE_EOF;
+    }
+    while (num_items_read != 3);
+
+    return 0;
+}
+
+
+int mm_read_mtx_array_size(FILE *f, int *M, int *N)
+{
+    char line[MM_MAX_LINE_LENGTH];
+    int num_items_read;
+    /* set return null parameter values, in case we exit with errors */
+    *M = *N = 0;
+
+    /* now continue scanning until you reach the end-of-comments */
+    do
+    {
+        if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
+            return MM_PREMATURE_EOF;
+    }while (line[0] == '%');
+
+    /* line[] is either blank or has M,N, nz */
+    if (sscanf(line, "%d %d", M, N) == 2)
+        return 0;
+
+    else /* we have a blank line */
+    do
+    {
+        num_items_read = fscanf(f, "%d %d", M, N);
+        if (num_items_read == EOF) return MM_PREMATURE_EOF;
+    }
+    while (num_items_read != 2);
+
+    return 0;
+}
+
+int mm_write_mtx_array_size(FILE *f, int M, int N)
+{
+    if (fprintf(f, "%d %d\n", M, N) != 2)
+        return MM_COULD_NOT_WRITE_FILE;
+    else
+        return 0;
+}
+
+
+
+/*-------------------------------------------------------------------------*/
+
+/******************************************************************/
+/* use when row[], col[], and val[]J, and val[] are already allocated */
+/******************************************************************/
+
+int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int row[], int col[],
+                         double val[], MM_typecode matcode)
+{
+    int i;
+    (void)M; (void)N;
+
+    if (mm_is_complex(matcode))
+    {
+        for (i=0; i<nz; i++)
+            if (fscanf(f, "%d %d %lg %lg", &row[i], &col[i], &val[2*i], &val[2*i+1])
+                != 4) return MM_PREMATURE_EOF;
+    }
+    else if (mm_is_real(matcode))
+    {
+        for (i=0; i<nz; i++)
+        {
+            if (fscanf(f, "%d %d %lg\n", &row[i], &col[i], &val[i])
+                != 3) return MM_PREMATURE_EOF;
+
+        }
+    }
+
+    else if (mm_is_pattern(matcode))
+    {
+        for (i=0; i<nz; i++)
+            if (fscanf(f, "%d %d", &row[i], &col[i])
+                != 2) return MM_PREMATURE_EOF;
+    }
+    else
+        return MM_UNSUPPORTED_TYPE;
+
+    return 0;
+
+}
+
+int mm_read_mtx_crd_entry(FILE *f, int *row, int *col,
+        double *real, double *imag, MM_typecode matcode)
+{
+    if (mm_is_complex(matcode))
+    {
+            if (fscanf(f, "%d %d %lg %lg", row, col, real, imag)
+                != 4) return MM_PREMATURE_EOF;
+    }
+    else if (mm_is_real(matcode))
+    {
+            if (fscanf(f, "%d %d %lg\n", row, col, real)
+                != 3) return MM_PREMATURE_EOF;
+
+    }
+
+    else if (mm_is_pattern(matcode))
+    {
+            if (fscanf(f, "%d %d", row, col) != 2) return MM_PREMATURE_EOF;
+    }
+    else
+        return MM_UNSUPPORTED_TYPE;
+
+    return 0;
+
+}
+
+
+/************************************************************************
+    mm_read_mtx_crd()  fills M, N, nz, array of values, and return
+                        type code, e.g. 'MCRS'
+
+                        if matrix is complex, values[] is of size 2*nz,
+                            (nz pairs of real/imaginary values)
+************************************************************************/
+
+int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **row, int **col,
+        double **val, MM_typecode *matcode)
+{
+    int ret_code;
+    FILE *f;
+
+    if (strcmp(fname, "stdin") == 0) f=stdin;
+    else
+    if ((f = fopen(fname, "r")) == NULL)
+        return MM_COULD_NOT_READ_FILE;
+
+
+    if ((ret_code = mm_read_banner(f, matcode)) != 0)
+        return ret_code;
+
+    if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) &&
+            mm_is_matrix(*matcode)))
+        return MM_UNSUPPORTED_TYPE;
+
+    if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
+        return ret_code;
+
+
+    *row = (int *)  malloc(*nz * sizeof(int));
+    *col = (int *)  malloc(*nz * sizeof(int));
+    *val = NULL;
+
+    if (mm_is_complex(*matcode))
+    {
+        *val = (double *) malloc(*nz * 2 * sizeof(double));
+        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *row, *col, *val,
+                *matcode);
+        if (ret_code != 0) return ret_code;
+    }
+    else if (mm_is_real(*matcode))
+    {
+        *val = (double *) malloc(*nz * sizeof(double));
+        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *row, *col, *val,
+                *matcode);
+        if (ret_code != 0) return ret_code;
+    }
+
+    else if (mm_is_pattern(*matcode))
+    {
+        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *row, *col, *val,
+                *matcode);
+        if (ret_code != 0) return ret_code;
+    }
+
+    if (f != stdin) fclose(f);
+    return 0;
+}
+
+int mm_write_banner(FILE *f, MM_typecode matcode)
+{
+    char *str = mm_typecode_to_str(matcode);
+    int ret_code;
+
+    ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
+    free(str);
+    if (ret_code !=2 )
+        return MM_COULD_NOT_WRITE_FILE;
+    else
+        return 0;
+}
+
+int mm_write_mtx_crd(char fname[], int M, int N, int nz, int row[], int col[],
+        double val[], MM_typecode matcode)
+{
+    FILE *f;
+    int i;
+
+    if (strcmp(fname, "stdout") == 0)
+        f = stdout;
+    else
+    if ((f = fopen(fname, "w")) == NULL)
+        return MM_COULD_NOT_WRITE_FILE;
+
+    /* print banner followed by typecode */
+    fprintf(f, "%s ", MatrixMarketBanner);
+    fprintf(f, "%s\n", mm_typecode_to_str(matcode));
+
+    /* print matrix sizes and nonzeros */
+    fprintf(f, "%d %d %d\n", M, N, nz);
+
+    /* print values */
+    if (mm_is_pattern(matcode))
+        for (i=0; i<nz; i++)
+            fprintf(f, "%d %d\n", row[i], col[i]);
+    else
+    if (mm_is_real(matcode))
+        for (i=0; i<nz; i++)
+            fprintf(f, "%d %d %20.16g\n", row[i], col[i], val[i]);
+    else
+    if (mm_is_complex(matcode))
+        for (i=0; i<nz; i++)
+            fprintf(f, "%d %d %20.16g %20.16g\n", row[i], col[i], val[2*i],
+                        val[2*i+1]);
+    else
+    {
+        if (f != stdout) fclose(f);
+        return MM_UNSUPPORTED_TYPE;
+    }
+
+    if (f !=stdout) fclose(f);
+
+    return 0;
+}
+
+
+/**
+*  Create a new copy of a string s.  mm_strdup() is a common routine, but
+*  not part of ANSI C, so it is included here.  Used by mm_typecode_to_str().
+*
+*/
+char *mm_strdup(const char *s)
+{
+        int len = strlen(s);
+        char *s2 = (char *) malloc((len+1)*sizeof(char));
+        return strcpy(s2, s);
+}
+
+char  *mm_typecode_to_str(MM_typecode matcode)
+{
+    char buffer[MM_MAX_LINE_LENGTH];
+    char *types[4];
+    char *mm_strdup(const char *);
+
+    /* check for MTX type */
+    if (mm_is_matrix(matcode))
+        types[0] = MM_MTX_STR;
+    else
+        return NULL;
+
+    /* check for CRD or ARR matrix */
+    if (mm_is_sparse(matcode))
+        types[1] = MM_SPARSE_STR;
+    else
+    if (mm_is_dense(matcode))
+        types[1] = MM_DENSE_STR;
+    else
+        return NULL;
+
+    /* check for element data type */
+    if (mm_is_real(matcode))
+        types[2] = MM_REAL_STR;
+    else
+    if (mm_is_complex(matcode))
+        types[2] = MM_COMPLEX_STR;
+    else
+    if (mm_is_pattern(matcode))
+        types[2] = MM_PATTERN_STR;
+    else
+    if (mm_is_integer(matcode))
+        types[2] = MM_INT_STR;
+    else
+        return NULL;
+
+
+    /* check for symmetry type */
+    if (mm_is_general(matcode))
+        types[3] = MM_GENERAL_STR;
+    else
+    if (mm_is_symmetric(matcode))
+        types[3] = MM_SYMM_STR;
+    else
+    if (mm_is_hermitian(matcode))
+        types[3] = MM_HERM_STR;
+    else
+    if (mm_is_skew(matcode))
+        types[3] = MM_SKEW_STR;
+    else
+        return NULL;
+
+    sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
+    return mm_strdup(buffer);
+
+}

drivers/mmio.h

+/*
+*   Matrix Market I/O library for ANSI C
+*
+*   See http://math.nist.gov/MatrixMarket for details.
+*
+*
+*/
+
+#ifndef MM_IO_H
+#define MM_IO_H
+
+#define MM_MAX_LINE_LENGTH 1025
+#define MatrixMarketBanner "%%MatrixMarket"
+#define MM_MAX_TOKEN_LENGTH 64
+
+typedef char MM_typecode[4];
+
+char *mm_typecode_to_str(MM_typecode matcode);
+
+int mm_read_banner(FILE *f, MM_typecode *matcode);
+int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz);
+int mm_read_mtx_array_size(FILE *f, int *M, int *N);
+
+int mm_write_banner(FILE *f, MM_typecode matcode);
+int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz);
+int mm_write_mtx_array_size(FILE *f, int M, int N);
+
+
+/********************* MM_typecode query fucntions ***************************/
+
+#define mm_is_matrix(typecode)     ((typecode)[0]=='M')
+
+#define mm_is_sparse(typecode)     ((typecode)[1]=='C')
+#define mm_is_coordinate(typecode) ((typecode)[1]=='C')
+#define mm_is_dense(typecode)      ((typecode)[1]=='A')
+#define mm_is_array(typecode)      ((typecode)[1]=='A')
+
+#define mm_is_complex(typecode)    ((typecode)[2]=='C')
+#define mm_is_real(typecode)       ((typecode)[2]=='R')
+#define mm_is_pattern(typecode)    ((typecode)[2]=='P')
+#define mm_is_integer(typecode)    ((typecode)[2]=='I')
+
+#define mm_is_symmetric(typecode)  ((typecode)[3]=='S')
+#define mm_is_general(typecode)    ((typecode)[3]=='G')
+#define mm_is_skew(typecode)       ((typecode)[3]=='K')
+#define mm_is_hermitian(typecode)  ((typecode)[3]=='H')
+
+int mm_is_valid(MM_typecode matcode);		/* too complex for a macro */
+
+
+/********************* MM_typecode modify fucntions ***************************/
+
+#define mm_set_matrix(typecode)      ((*typecode)[0]='M')
+#define mm_set_coordinate(typecode)  ((*typecode)[1]='C')
+#define mm_set_array(typecode)       ((*typecode)[1]='A')
+#define mm_set_dense(typecode)       mm_set_array(typecode)
+#define mm_set_sparse(typecode)      mm_set_coordinate(typecode)
+
+#define mm_set_complex(typecode)     ((*typecode)[2]='C')
+#define mm_set_real(typecode)        ((*typecode)[2]='R')
+#define mm_set_pattern(typecode)     ((*typecode)[2]='P')
+#define mm_set_integer(typecode)     ((*typecode)[2]='I')
+
+
+#define mm_set_symmetric(typecode)  ((*typecode)[3]='S')
+#define mm_set_general(typecode)    ((*typecode)[3]='G')
+#define mm_set_skew(typecode)       ((*typecode)[3]='K')
+#define mm_set_hermitian(typecode)  ((*typecode)[3]='H')
+
+#define mm_clear_typecode(typecode) ((*typecode)[0]=(*typecode)[1]= \
+                                     (*typecode)[2]=' ',(*typecode)[3]='G')
+
+#define mm_initialize_typecode(typecode) mm_clear_typecode(typecode)
+
+
+/********************* Matrix Market error codes ***************************/
+
+
+#define MM_COULD_NOT_READ_FILE	11
+#define MM_PREMATURE_EOF	12
+#define MM_NOT_MTX		13
+#define MM_NO_HEADER		14
+#define MM_UNSUPPORTED_TYPE	15
+#define MM_LINE_TOO_LONG	16
+#define MM_COULD_NOT_WRITE_FILE	17
+
+
+/******************** Matrix Market internal definitions ********************
+
+   MM_matrix_typecode: 4-character sequence
+        ojbect          sparse/   data        storage
+        dense           type      scheme
+
+   string position:      [0]        [1]			[2]         [3]
+
+   Matrix typecode:  M(atrix)   C(oord)     R(eal)      G(eneral)
+                     A(array)	C(omplex)   H(ermitian)
+                     P(attern)  S(ymmetric)
+                     I(nteger)	K(kew)
+
+ ***********************************************************************/
+
+#define MM_MTX_STR        "matrix"
+#define MM_ARRAY_STR      "array"
+#define MM_DENSE_STR      "array"
+#define MM_COORDINATE_STR "coordinate"
+#define MM_SPARSE_STR     "coordinate"
+#define MM_COMPLEX_STR    "complex"
+#define MM_REAL_STR       "real"
+#define MM_INT_STR        "integer"
+#define MM_GENERAL_STR    "general"
+#define MM_SYMM_STR       "symmetric"
+#define MM_HERM_STR       "hermitian"
+#define MM_SKEW_STR       "skew-symmetric"
+#define MM_PATTERN_STR    "pattern"
+
+
+/*  high level routines */
+int mm_write_mtx_crd(char fname[], int M, int N, int nz, int Itab[], int Jtab[],
+                     double val[], MM_typecode matcode);
+int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int Itab[], int Jtab[],
+                         double val[], MM_typecode matcode);
+int mm_read_mtx_crd_entry(FILE *f, int *Itab, int *Jtab, double *real, double *img,
+                          MM_typecode matcode);
+
+int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
+                double **val_, int **I_, int **J_);
+
+
+
+#endif

drivers/readhb.c

+/**
+ * @file readhb.c
+ *
+ *  $COPYRIGHTS$
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Pierre Ramet
+ * @author Xavier Lacoste
+ * @date 2011-11-11
+ *
+ **/
+#include <stdio.h>
+#include "common.h"
+#include "spm_drivers.h"
+#include "drivers/iohb.h"
+
+/**
+ * ******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * readHB - Interface to the Harwell-Boeing C driver (iohb.c)
+ *
+ *******************************************************************************
+ *
+ * @param[in] filename
+ *          The file containing the matrix.
+ *
+ * @param[in] csc
+ *          At exit, contains the matrix in csc format.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been read successfully
+ *      \retval PASTIX_ERR_IO if a problem occured in the Harwell Boeing driver
+ *      \retval PASTIX_ERR_BADPARAMETER if the matrix is no in a supported format
+ *
+ *******************************************************************************/
+int
+readHB( const char   *filename,
+        pastix_csc_t *csc )
+{
+    int M, N, nz, nrhs;
+
+    /* Harwell Boeing is a variant of RSA */
+    csc->fmttype = PastixCSC;
+    csc->dof     = 1;
+    csc->loc2glob= NULL;
+
+    /* Read header informations */
+    {
+        char *Type;
+        Type = malloc(4*sizeof(char));
+        Type[0] ='a';
+
+        readHB_info(filename, &M, &N, &nz, &Type, &nrhs);
+
+        if ( M != N ) {
+            fprintf(stderr, "readHB: PaStiX does not support non square matrices (m=%d, N=%d\n", M, N);
+            return PASTIX_ERR_BADPARAMETER;
+        }
+
+        csc->gN   = M;
+        csc->n    = M;
+        csc->gnnz = nz;
+        csc->nnz  = nz;
+
+        /* Check float type */
+        switch( Type[0] ) {
+        case 'C':
+        case 'c':
+            csc->flttype = PastixComplex64;
+            break;
+        case 'R':
+        case 'r':
+            csc->flttype = PastixDouble;
+            break;
+        case 'P':
+        case 'p':
+            csc->flttype = PastixPattern;
+            break;
+        default:
+            fprintf(stderr, "readhb: Floating type unknown (%c)\n", Type[0]);
+            return PASTIX_ERR_BADPARAMETER;
+        }
+
+        /* Check Symmetry */
+        switch( Type[1] ) {
+        case 'S':
+        case 's':
+            csc->mtxtype = PastixSymmetric;
+            break;
+        case 'H':
+        case 'h':
+            csc->mtxtype = PastixHermitian;
+            assert( csc->flttype == PastixDouble );
+            break;
+        case 'U':
+        case 'u':
+        default:
+            csc->mtxtype = PastixGeneral;
+        }
+        free(Type);
+    }
+
+    /* Read the matrix and its values */
+    {
+        int    *colptr, *rowind;
+        int     rc;
+
+        rc = readHB_newmat_double( filename, &M, &N, &nz,
+                                   &colptr, &rowind, (double**)(&(csc->values)) );
+
+        if (rc == 0) {
+            fprintf(stderr, "readhb: Error in reading the HB matrix values\n");
+            return PASTIX_ERR_IO;
+        }
+
+        /* Move the colptr/rowind from int to pastix_int_t if different sizes */
+        csc->colptr = spmIntConvert(csc->n+1, colptr);
+        csc->rowptr = spmIntConvert(csc->nnz, rowind);
+    }
+    return PASTIX_SUCCESS;
+}

drivers/readijv.c

+/**
+ * @file readijv.c
+ *
+ *  $COPYRIGHTS$
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Pierre Ramet
+ * @author Xavier Lacoste
+ * @date 2011-11-11
+ *
+ **/
+#include <stdio.h>
+#include <stdlib.h>
+#include "common.h"
+#include "spm_drivers.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * threeFilesReadHeader - Read header from three file IJV format.
+ *
+ *******************************************************************************
+ *
+ * @param[in] infile
+ *          The opened header file
+ *
+ * @param[out] Nrow
+ *          At exit, contains the number of rows of the matrix.
+ *
+ * @param[out] Ncol
+ *          At exit, contains the number of columns of the matrix.
+ *
+ * @param[out] Nnzero
+ *          At exit, contains the number of non zero entries of the matrix.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the information has been read successfully
+ *      \retval PASTIX_ERR_BADPARAMETER if the header has a wrong format
+ *
+ *******************************************************************************/
+int
+threeFilesReadHeader(FILE         *infile,
+                     pastix_int_t *Nrow,
+                     pastix_int_t *Ncol,
+                     pastix_int_t *Nnzero)
+{
+    long temp1,temp2,temp3;
+
+    /* ncol nrow nnzero */
+    if (fscanf(infile, "%ld %ld %ld\n", &temp1, &temp2, &temp3) != 3) {
+        Nrow = Ncol = Nnzero = 0;
+        fprintf(stderr, "readijv: Wrong format in header file\n");
+        return PASTIX_ERR_BADPARAMETER;
+    }
+    *Nrow   = (pastix_int_t)temp1;
+    *Ncol   = (pastix_int_t)temp2;
+    *Nnzero = (pastix_int_t)temp3;
+
+    return PASTIX_SUCCESS;
+}
+
+/**
+ * ******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * readIJV - Read matrix from three files IJV
+ *
+ * header file is "filename"/header
+ * columns file is "filename"/ia_threeFiles
+ * rows file is "filename"/ja_threeFiles
+ * values file is "filename"/ra_threeFiles
+ *
+ *******************************************************************************
+ *
+ * @param[in] dirname
+ *          Directory that contains the files.
+ *
+ * @param[out] csc
+ *          At exit, contains the matrix in ijv format.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been read successfully
+ *      \retval PASTIX_ERR_IO if a problem occurs while reading the files
+ *      \retval PASTIX_ERR_BADPARAMETER if a problem occurs while opening the
+ *              files
+ *
+ *******************************************************************************/
+int
+readIJV( const char   *dirname,
+         pastix_csc_t *csc )
+{
+
+    FILE *iafile, *jafile, *rafile;
+    FILE *hdrfile;
+    char *filename;
+    pastix_int_t *tempcol;
+    pastix_int_t *temprow;
+    double       *tempval;
+    pastix_int_t  i, Nrow, Ncol, Nnzero;
+
+    filename = malloc(strlen(dirname)+10);
+
+    csc->flttype = PastixDouble;
+    csc->mtxtype = PastixGeneral;
+    csc->fmttype = PastixIJV;
+    csc->dof     = 1;
+    csc->loc2glob= NULL;
+
+    /* Read the header information */
+    {
+        sprintf(filename,"%s/header",dirname);
+        hdrfile = fopen (filename,"r");
+        if (hdrfile == NULL)
+        {
+            fprintf(stderr,"readijv: Cannot open the header file (%s)\n", filename);
+            return PASTIX_ERR_BADPARAMETER;
+        }
+        threeFilesReadHeader(hdrfile, &Nrow, &Ncol, &Nnzero);
+        fclose(hdrfile);
+    }
+
+    csc->gN      = Ncol;
+    csc->n       = Ncol;
+    csc->gnnz    = Nnzero;
+    csc->nnz     = Nnzero;
+    csc->colptr = (pastix_int_t *) malloc(Nnzero*sizeof(pastix_int_t));
+    csc->rowptr = (pastix_int_t *) malloc(Nnzero*sizeof(pastix_int_t));
+    csc->values = (double *)       malloc(Nnzero*sizeof(double));
+
+    /* Open the 3 files */
+    sprintf(filename,"%s/ia_threeFiles",dirname);
+    iafile = fopen(filename,"r");
+    if (iafile == NULL)
+    {
+        fprintf(stderr,"readijv: Cannot open the ia file (%s)\n", filename);
+        return PASTIX_ERR_BADPARAMETER;
+    }
+
+    sprintf(filename,"%s/ja_threeFiles",dirname);
+    jafile = fopen(filename,"r");
+    if (jafile == NULL)
+    {
+        fprintf(stderr,"readijv: Cannot open the ja file (%s)\n", filename);
+        fclose(iafile);
+        return PASTIX_ERR_BADPARAMETER;
+    }
+
+    sprintf(filename,"%s/ra_threeFiles",dirname);
+    rafile = fopen(filename,"r");
+    if (rafile == NULL)
+    {
+        fprintf(stderr,"readijv: Cannot open the ra file (%s)\n", filename);
+        fclose(iafile);
+        fclose(jafile);
+        return PASTIX_ERR_BADPARAMETER;
+    }
+
+    /* Read the files */
+    tempcol = csc->colptr;
+    temprow = csc->rowptr;
+    tempval = csc->values;
+
+    for (i=0; i<Nnzero; i++, tempcol++, temprow++, tempval++)
+    {
+        long temp1, temp2;
+        double temp3;
+
+        if (( 1 != fscanf(iafile,"%ld\n", &temp1)) ||
+            ( 1 != fscanf(jafile,"%ld\n", &temp2)) ||
+            ( 1 != fscanf(rafile,"%le\n", &temp3)) )
+        {
+            fprintf(stderr, "ERROR: reading matrix\n");
+            return PASTIX_ERR_IO;
+        }
+        *temprow = (pastix_int_t)temp1;
+        *tempcol = (pastix_int_t)temp2;
+        *tempval = temp3;
+    }
+    fclose(iafile);
+    fclose(jafile);
+    fclose(rafile);
+
+    return PASTIX_SUCCESS;
+}

drivers/readmm.c

+/**
+ * @file readmm.c
+ *
+ *  $COPYRIGHTS$
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Pierre Ramet
+ * @author Xavier Lacoste
+ * @author Theophile Terraz
+ * @date 2011-11-11
+ *
+ **/
+#include "common.h"
+#include "spm_drivers.h"
+#include "drivers/mmio.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * z_readMM - Read the data part of a complex matrix in Matrix Market file.
+ * For more information about matrix market format see mmio.c/mmio.h
+ *
+ *******************************************************************************
+ *
+ * @param[in] file
+ *          The file opened in readMM which contains the matrix stored in Matrix
+ *          Market format.
+ *
+ * @param[in,out] csc
+ *          At exit, the data of the matrix are stored in the csc structure.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been read successfully
+ *      \retval PASTIX_ERR_IO if a problem occured in the RSA driver
+ *
+ *******************************************************************************/
+int
+z_readMM( FILE *file,
+          pastix_csc_t *csc )
+{
+    pastix_complex64_t *valptr;
+    pastix_int_t *colptr;
+    pastix_int_t *rowptr;
+    pastix_int_t i;
+    long row, col;
+    double re, im;
+
+    csc->values = malloc( csc->nnz * sizeof(pastix_complex64_t) );
+
+    colptr = csc->colptr;
+    rowptr = csc->rowptr;
+    valptr = (pastix_complex64_t*)(csc->values);
+
+    for (i=0; i<csc->nnz; i++, colptr++, rowptr++, valptr++)
+    {
+        if (4 != fscanf(file,"%ld %ld %lg %lg\n", &row, &col, &re, &im))
+        {
+            fprintf(stderr, "readmm: erro while reading matrix file (line %ld)\n", (long)i);
+            return PASTIX_ERR_IO;
+        }
+
+        *rowptr = (pastix_int_t)row;
+        *colptr = (pastix_int_t)col;
+        *valptr = (pastix_complex64_t)(re + im * I);
+    }
+
+    return PASTIX_SUCCESS;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * d_readMM - Read the data part of a real matrix in Matrix Market file.
+ * For more information about matrix market format see mmio.c/mmio.h
+ *
+ *******************************************************************************
+ *
+ * @param[in] file
+ *          The file opened in readMM which contains the matrix stored in Matrix
+ *          Market format.
+ *
+ * @param[in,out] csc
+ *          At exit, the data of the matrix are stored in the csc structure.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been read successfully
+ *      \retval PASTIX_ERR_IO if a problem occured in the RSA driver
+ *
+ *******************************************************************************/
+int
+d_readMM( FILE *file,
+          pastix_csc_t *csc )
+{
+    double       *valptr;
+    pastix_int_t *colptr;
+    pastix_int_t *rowptr;
+    pastix_int_t i;
+    long row, col;
+    double re;
+
+    csc->values = malloc( csc->nnz * sizeof(double) );
+
+    colptr = csc->colptr;
+    rowptr = csc->rowptr;
+    valptr = (double*)(csc->values);
+
+    for (i=0; i<csc->nnz; i++, colptr++, rowptr++, valptr++)
+    {
+        if (3 != fscanf(file,"%ld %ld %lg\n", &row, &col, &re))
+        {
+            fprintf(stderr, "readmm: erro while reading matrix file (line %ld)\n", (long)i);
+            return PASTIX_ERR_IO;
+        }
+
+        *rowptr = (pastix_int_t)row;
+        *colptr = (pastix_int_t)col;
+        *valptr = re;
+    }
+
+    return PASTIX_SUCCESS;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * p_readMM - Read the data part of a pattern matrix in Matrix Market file.
+ * For more information about matrix market format see mmio.c/mmio.h
+ *
+ *******************************************************************************
+ *
+ * @param[in] file
+ *          The file opened in readMM which contains the matrix stored in Matrix
+ *          Market format.
+ *
+ * @param[in,out] csc
+ *          At exit, the data of the matrix are stored in the csc structure.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been read successfully
+ *      \retval PASTIX_ERR_IO if a problem occured in the RSA driver
+ *
+ *******************************************************************************/
+int
+p_readMM( FILE *file,
+          pastix_csc_t *csc )
+{
+    pastix_int_t *colptr;
+    pastix_int_t *rowptr;
+    pastix_int_t i;
+    long row, col;
+
+    csc->values = NULL;
+
+    colptr = csc->colptr;
+    rowptr = csc->rowptr;
+
+    for (i=0; i<csc->nnz; i++, colptr++, rowptr++)
+    {
+        if (2 != fscanf(file,"%ld %ld\n", &row, &col))
+        {
+            fprintf(stderr, "readmm: erro while reading matrix file (line %ld)\n", (long)i);
+            return PASTIX_ERR_IO;
+        }
+
+        *rowptr = (pastix_int_t)row;
+        *colptr = (pastix_int_t)col;
+    }
+
+    return PASTIX_SUCCESS;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * readMM - Read a matrix in Matrix Market fill. This corresponds to
+ * IJV format with (%d %d[ %lf[ %lf]]) format per line.
+ * For more information about matrix market format see mmio.c/mmio.h
+ *
+ *******************************************************************************
+ *
+ * @param[in] filename
+ *          The filename that contains the matrix stored in Matrix Market format.
+ *
+ * @param[in] csc
+ *          At exit, contains the matrix in csc format.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been read successfully
+ *      \retval PASTIX_ERR_IO if a problem occured in the RSA driver
+ *      \retval PASTIX_ERR_BADPARAMETER if the matrix is no in a supported format
+ *
+ *******************************************************************************/
+int
+readMM( const char   *filename,
+        pastix_csc_t *csc )
+{
+    MM_typecode matcode;
+    FILE *file;
+    int rc;
+
+    file = fopen(filename,"r");
+    if (file == NULL)
+    {
+        fprintf(stderr,"readmm: Cannot open the file (%s)\n", filename);
+        return PASTIX_ERR_BADPARAMETER;
+    }
+
+    if (mm_read_banner(file, &matcode) != 0)
+    {
+        fprintf(stderr,"readmm: Could not process Matrix Market banner.\n");
+        return PASTIX_ERR_IO;
+    }
+
+    /* Float values type */
+    if (mm_is_complex(matcode)) {
+        csc->flttype = PastixComplex64;
+    }
+    else if (mm_is_real(matcode)) {
+        csc->flttype = PastixDouble;
+    }
+    else if (mm_is_pattern(matcode)) {
+        csc->flttype = PastixPattern;
+    }
+    else {
+        fprintf(stderr,"readmm: Unsupported type of matrix.\n");
+        return PASTIX_ERR_BADPARAMETER;
+    }
+
+    /* Matrix structure */
+    if (mm_is_general(matcode)) {
+        csc->mtxtype = PastixGeneral;
+    }
+    else if (mm_is_symmetric(matcode)) {
+        csc->mtxtype = PastixSymmetric;
+    }
+    else if (mm_is_hermitian(matcode)) {
+        csc->mtxtype = PastixHermitian;
+    }
+    else {
+        fprintf(stderr,"readmm: Unsupported type of matrix.\n");
+        return PASTIX_ERR_BADPARAMETER;
+    }
+
+    csc->fmttype = PastixIJV;
+    csc->dof     = 1;
+    csc->loc2glob= NULL;
+
+    /* Read the size */
+    {
+        int m, n, nnz;
+        if (mm_read_mtx_crd_size(file, &m, &n, &nnz) != 0) {
+            fprintf(stderr, "readmm: error while reading matrix sizes\n");
+            return PASTIX_ERR_IO;
+        }
+
+        csc->gN   = n;
+        csc->n    = n;
+        csc->gnnz = nnz;
+        csc->nnz  = nnz;
+    }
+
+    csc->colptr = (pastix_int_t*)malloc(csc->nnz * sizeof(pastix_int_t));
+    csc->rowptr = (pastix_int_t*)malloc(csc->nnz * sizeof(pastix_int_t));
+
+    switch( csc->flttype ) {
+    case PastixComplex64:
+        rc = z_readMM(file,csc);
+        break;
+
+    case PastixDouble:
+        rc = d_readMM(file,csc);
+        break;
+
+    case PastixPattern:
+    default:
+        rc = p_readMM(file,csc);
+    }
+
+    fclose(file);
+    return rc;
+}

drivers/readrsa.c

+/**
+ * @file readrsa.c
+ *
+ *  $COPYRIGHTS$
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Pierre Ramet
+ * @author Xavier Lacoste
+ * @date 2011-11-11
+ *
+ **/
+#include "common.h"
+#include "spm_drivers.h"
+
+/**
+ * Function: FORTRAN_CALL(wreadmtc)
+ *
+ * Declaration of the wreadmtc fortran function
+ * defined in skitf.f.
+ *
+ * Parameters:
+ *    tmp1      - Maximum number of column                                    (INPUT)
+ *    tmp2      - Maximum number of non zeros                                 (INPUT)
+ *    tmp3      - job to be done (see skitf file)                             (INPUT)
+ *    filename  - Path to the file to read from                               (INPUT)
+ *    len       - length of *filname*                                         (INPUT)
+ *    val       - Values of the elements of the matrix.                       (OUTPUT)
+ *    row       - Rows of the elements of the matrix.                         (OUTPUT)
+ *    col       - Index of first element of each column in *row* and *val*    (OUTPUT)
+ *    crhs      - Right hand side(s).                                         (OUTPUT)
+ *    nrhs      - Number of right hand side(s).                               (OUTPUT)
+ *    RhsType   - Right hand side type                                        (OUTPUT)
+ *    tmpNrow   - Number of rows.                                             (OUTPUT)
+ *    tmpNcol   - Number of columns.                                          (OUTPUT)
+ *    tmpNnzero - Number of non zeros.                                        (OUTPUT)
+ *    title     - name of the matrix.                                         (OUTPUT)
+ *    key       - key of the matrix (see skitf.f)                             (OUTPUT)
+ *    Type      - Type of the matrix                                          (OUTPUT)
+ *    ierr      - Error return value                                          (OUTPUT)
+ *
+ */
+void
+FC_GLOBAL(wreadmtc,WREADMTC)(int        *tmp1,
+                             int        *tmp2,
+                             int        *tmp3,
+                             const char *filename,
+                             int        *len,
+                             double     *val,
+                             int        *row,
+                             int        *col,
+                             double     *crhs,
+                             int        *nrhs,
+                             char       *RhsType,
+                             int        *tmpNrow,
+                             int        *tmpNcol,
+                             int        *tmpNnzero,
+                             char       *title,
+                             char       *key,
+                             char       *Type,
+                             int        *ierr);
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * readRSAHeader - Read the header structure of a RSA file
+ *
+ *******************************************************************************
+ *
+ * @param[in] filename
+ *          The filename that contains the matrix stored in RSA format. This is
+ *          an interface to the wreadmtx fortran fucntion provided within
+ *          SparseKit.
+ *
+ * @param[out] N
+ *          At exit, contains the number of rows/columns of the matrix.
+ *
+ * @param[out] Nnzero
+ *          At exit, contains the number of non zero entries of the matrix.
+ *
+ * @param[out] Type
+ *          At exit, contains the type of the matrix.
+ *
+ * @param[out] RhsType
+ *          At exit, contains the type of the right hand side.
+ *
+ *******************************************************************************/
+void
+readRSAHeader( const char *filename,
+               int        *N,
+               int        *Nnz,
+               char       *Type,
+               char       *RhsType )
+{
+    int     tmp;
+    int    *col = NULL;
+    int    *row = NULL;
+    char    title[72+1];
+    char    key[8+1];
+    int     nrhs;
+    int     len;
+    int     ierr;
+    double *val  = NULL;
+    double *crhs = NULL;
+    int     M;
+
+    len = strlen(filename);
+    tmp = 0;
+
+    FC_GLOBAL(wreadmtc,WREADMTC)
+        (&tmp, &tmp, &tmp, filename, &len, val, row, col, crhs, &nrhs,
+         RhsType, &M, N, Nnz, title, key, Type, &ierr );
+
+    if(ierr != 0) {
+        fprintf(stderr, "cannot read matrix (job=0)\n");
+    }
+
+    if ( M != (*N) )
+    {
+        fprintf(stderr,"ERROR : (M != N)\n");
+        exit(EXIT_FAILURE);
+    }
+
+    Type[3] = '\0';
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup pastix_csc_driver
+ *
+ * readRSA - Read a RSA matrix file. This driver reads only real matrices, and
+ * does not support complex matrices.
+ * The matrix is returned in double, convert it to real if needed through TODO
+ *
+ *******************************************************************************
+ *
+ * @param[in] filename
+ *          The filename that contains the matrix stored in RSA format. This is
+ *          an interface to the wreadmtx fortran function provided within
+ *          SparseKit.
+ *
+ * @param[in] csc
+ *          At exit, contains the matrix in csc format.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *      \retval PASTIX_SUCCESS if the matrix has been read successfully
+ *      \retval PASTIX_ERR_IO if a problem occured in the RSA driver
+ *      \retval PASTIX_ERR_BADPARAMETER if the matrix is no in a supported format
+ *
+ *******************************************************************************/
+int
+readRSA( const char   *filename,
+         pastix_csc_t *csc )
+{
+    char    Type[4];
+    char    RhsType[4];
+    int     tmp;
+    char    title[72+1];
+    char    key[8+1];
+    int     nrhs;
+    int     len;
+    int     ierr;
+    double *crhs=NULL;
+    int     M, N, Nnz;
+    int    *tmpcolptr;
+    int    *tmprows;
+
+    readRSAHeader(filename, &N, &Nnz, Type, RhsType );
+
+    switch( Type[1] ){
+    case 'S':
+    case 's':
+        csc->mtxtype = PastixSymmetric;
+        break;
+    case 'H':
+    case 'h':
+        csc->mtxtype = PastixHermitian;
+        /**
+         * We should not arrive here, since the fortran driver is not able to
+         * read complex matrices
+         */
+        fprintf(stderr,"readrsa: Unsupported Complex.\n");
+        return PASTIX_ERR_BADPARAMETER;
+    case 'U':
+    case 'u':
+        csc->mtxtype = PastixGeneral;
+        break;
+    default:
+        fprintf(stderr,"readrsa: Unsupported type of matrix.\n");
+        return PASTIX_ERR_BADPARAMETER;
+    }
+
+    csc->flttype = PastixDouble;
+    csc->fmttype = PastixCSC;
+    csc->gN      = N;
+    csc->n       = N;
+    csc->gnnz    = Nnz;
+    csc->nnz     = Nnz;
+    csc->dof     = 1;
+    csc->loc2glob= NULL;
+
+    tmpcolptr = (int*) malloc( (N+1) * sizeof(int) );
+    assert( tmpcolptr );
+
+    tmprows = (int*) malloc( Nnz * sizeof(int) );
+    assert( tmprows );
+
+    /* RSA reads only double */
+    csc->values = (double*) malloc( Nnz * sizeof(double) );
+    assert( csc->values );
+
+    len  = strlen(filename);
+    tmp  = 2;
+    nrhs = 0;
+
+    FC_GLOBAL(wreadmtc,WREADMTC)
+        (&N, &Nnz, &tmp, filename, &len, csc->values, tmprows, tmpcolptr, crhs,
+         &nrhs, RhsType, &M, &N, &Nnz, title, key, Type, &ierr );
+
+    assert( (tmpcolptr[N]-tmpcolptr[0]) == Nnz );
+
+    csc->colptr  = spmIntConvert( N+1, tmpcolptr );
+    csc->rowptr  = spmIntConvert( Nnz, tmprows );
+
+    RhsType[0] = '\0';
+    if(ierr != 0) {
+        fprintf(stderr, "cannot read matrix (job=2)\n");
+        return PASTIX_ERR_IO;
+    }
+
+    return PASTIX_SUCCESS;
+}

integer_sort.c

+/* This file is part of the Scotch distribution. It does
+** not have the stardard Scotch header with the INRIA
+** copyright notice because it is a very slight adaptation
+** of the qsort routine of glibc 2.4, taylored to match
+** Scotch needs. As Scotch is distributed according to the
+** CeCILL-C license, which is LGPL-compatible, no further
+** notices are required.
+*/
+
+/* Copyright (C) 1991,1992,1996,1997,1999,2004 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Written by Douglas C. Schmidt (schmidt@ics.uci.edu).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+/* If you consider tuning this algorithm, you should consult first:
+   Engineering a sort function; Jon Bentley and M. Douglas McIlroy;
+   Software - Practice and Experience; Vol. 23 (11), 1249-1265, 1993.  */
+
+
+/**
+ * @file integer_sort.c
+ *
+ * File to include to create integer sort function using qsort based
+ * algorithm. DO NOT compile directly.
+ *
+ * @author François Pellegrini
+ *
+*/
+
+#ifndef MAX_THRESH
+#define MAX_THRESH 6
+
+#define max_thresh                  (MAX_THRESH * INTSORTSIZE) /* Variable turned into constant */
+
+/* Stack node declarations used to store unfulfilled partition obligations. */
+typedef struct
+  {
+    char *lo;
+    char *hi;
+  } stack_node;
+
+/* The next 4 #defines implement a very fast in-line stack abstraction. */
+/* The stack needs log (total_elements) entries (we could even subtract
+   log(MAX_THRESH)).  Since total_elements has type size_t, we get as
+   upper bound for log (total_elements):
+   bits per byte (CHAR_BIT) * sizeof(size_t).  */
+#define STACK_SIZE	(CHAR_BIT * sizeof (pastix_int_t))
+#define PUSH(low, high)	((void) ((top->lo = (low)), (top->hi = (high)), ++top))
+#define	POP(low, high)	((void) (--top, (low = top->lo), (high = top->hi)))
+#define	STACK_NOT_EMPTY	(stack < top)
+
+#endif /* MAX_THRESH */
+
+/* Order size using quicksort.  This implementation incorporates
+   four optimizations discussed in Sedgewick:
+
+   1. Non-recursive, using an explicit stack of pointer that store the
+      next array partition to sort.  To save time, this maximum amount
+      of space required to store an array of SIZE_MAX is allocated on the
+      stack.  Assuming a 32-bit (64 bit) integer for size_t, this needs
+      only 32 * sizeof(stack_node) == 256 bytes (for 64 bit: 1024 bytes).
+      Pretty cheap, actually.
+
+   2. Chose the pivot element using a median-of-three decision tree.
+      This reduces the probability of selecting a bad pivot value and
+      eliminates certain extraneous comparisons.
+
+   3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving
+      insertion sort to order the MAX_THRESH items within each partition.
+      This is a big win, since insertion sort is faster for small, mostly
+      sorted array segments.
+
+   4. The larger of the two sub-partitions is always pushed onto the
+      stack first, with the algorithm then concentrating on the
+      smaller partition.  This *guarantees* no more than log (total_elems)
+      stack size is needed (actually O(1) in this case)!  */
+
+/* To be defined :
+** INTSORTNAME : Name of function
+** INTSORTSIZE : Size of elements to sort
+** INTSORTSWAP : Swapping macro
+** INTSORTCMP  : Comparison function
+*/
+
+void
+INTSORTNAME (
+void * const                pbase,                /*+ Array to sort             +*/
+const pastix_int_t          total_elems)          /*+ Number of entries to sort +*/
+{
+  register char *base_ptr = (char *) pbase;
+
+  if (total_elems == 0)
+    /* Avoid lossage with unsigned arithmetic below.  */
+    return;
+
+  if (total_elems > MAX_THRESH)
+    {
+      char *lo = base_ptr;
+      char *hi = &lo[INTSORTSIZE * (total_elems - 1)];
+      stack_node stack[STACK_SIZE];
+      stack_node *top = stack;
+
+      PUSH (NULL, NULL);
+
+      while (STACK_NOT_EMPTY)
+        {
+          char *left_ptr;
+          char *right_ptr;
+
+          /* Select median value from among LO, MID, and HI. Rearrange
+             LO and HI so the three values are sorted. This lowers the
+             probability of picking a pathological pivot value and
+             skips a comparison for both the LEFT_PTR and RIGHT_PTR in
+             the while loops. */
+
+          char *mid = lo + INTSORTSIZE * ((hi - lo) / INTSORTSIZE >> 1);
+
+          if (INTSORTCMP ((void *) mid, (void *) lo))
+            INTSORTSWAP (mid, lo);
+          if (INTSORTCMP ((void *) hi, (void *) mid))
+            INTSORTSWAP (mid, hi);
+          else
+            goto jump_over;
+          if (INTSORTCMP ((void *) mid, (void *) lo))
+            INTSORTSWAP (mid, lo);
+        jump_over:;
+
+          left_ptr  = lo + INTSORTSIZE;
+          right_ptr = hi - INTSORTSIZE;
+
+          /* Here's the famous ``collapse the walls'' section of quicksort.
+             Gotta like those tight inner loops!  They are the main reason
+             that this algorithm runs much faster than others. */
+          do
+            {
+              while (INTSORTCMP ((void *) left_ptr, (void *) mid))
+                left_ptr += INTSORTSIZE;
+
+              while (INTSORTCMP ((void *) mid, (void *) right_ptr))
+                right_ptr -= INTSORTSIZE;
+
+              if (left_ptr < right_ptr)
+                {
+                  INTSORTSWAP (left_ptr, right_ptr);
+                  if (mid == left_ptr)
+                    mid = right_ptr;
+                  else if (mid == right_ptr)
+                    mid = left_ptr;
+                  left_ptr += INTSORTSIZE;
+                  right_ptr -= INTSORTSIZE;
+                }
+              else if (left_ptr == right_ptr)
+                {
+                  left_ptr += INTSORTSIZE;
+                  right_ptr -= INTSORTSIZE;
+                  break;
+                }
+            }
+          while (left_ptr <= right_ptr);
+
+          /* Set up pointers for next iteration.  First determine whether
+             left and right partitions are below the threshold size.  If so,
+             ignore one or both.  Otherwise, push the larger partition's
+             bounds on the stack and continue sorting the smaller one. */
+
+          if ((size_t) (right_ptr - lo) <= max_thresh)
+            {
+              if ((size_t) (hi - left_ptr) <= max_thresh)
+                /* Ignore both small partitions. */
+                POP (lo, hi);
+              else
+                /* Ignore small left partition. */
+                lo = left_ptr;
+            }
+          else if ((size_t) (hi - left_ptr) <= max_thresh)
+            /* Ignore small right partition. */
+            hi = right_ptr;
+          else if ((right_ptr - lo) > (hi - left_ptr))
+            {
+              /* Push larger left partition indices. */
+              PUSH (lo, right_ptr);
+              lo = left_ptr;
+            }
+          else
+            {
+              /* Push larger right partition indices. */
+              PUSH (left_ptr, hi);
+              hi = right_ptr;
+            }
+        }
+    }
+
+  /* Once the BASE_PTR array is partially sorted by quicksort the rest
+     is completely sorted using insertion sort, since this is efficient
+     for partitions below MAX_THRESH size. BASE_PTR points to the beginning
+     of the array to sort, and END_PTR points at the very last element in
+     the array (*not* one beyond it!). */
+
+#define min(x, y) ((x) < (y) ? (x) : (y))
+
+  {
+    char *const end_ptr = &base_ptr[INTSORTSIZE * (total_elems - 1)];
+    char *tmp_ptr = base_ptr;
+    char *thresh = min(end_ptr, base_ptr + max_thresh);
+    register char *run_ptr;
+
+    /* Find smallest element in first threshold and place it at the
+       array's beginning.  This is the smallest array element,
+       and the operation speeds up insertion sort's inner loop. */
+
+    for (run_ptr = tmp_ptr + INTSORTSIZE; run_ptr <= thresh; run_ptr += INTSORTSIZE)
+      if (INTSORTCMP ((void *) run_ptr, (void *) tmp_ptr))
+        tmp_ptr = run_ptr;
+
+    if (tmp_ptr != base_ptr)
+      INTSORTSWAP (tmp_ptr, base_ptr);
+
+    /* Insertion sort, running from left-hand-side up to right-hand-side.  */
+
+    run_ptr = base_ptr + INTSORTSIZE;
+    while ((run_ptr += INTSORTSIZE) <= end_ptr)
+      {
+        tmp_ptr = run_ptr - INTSORTSIZE;
+        while (INTSORTCMP ((void *) run_ptr, (void *) tmp_ptr))
+          tmp_ptr -= INTSORTSIZE;
+
+        tmp_ptr += INTSORTSIZE;
+        if (tmp_ptr != run_ptr)
+          {
+            char *trav;
+
+            trav = run_ptr + INTSORTSIZE;
+            while (--trav >= run_ptr)
+              {
+                char c = *trav;
+                char *hi, *lo;
+
+                for (hi = lo = trav; (lo -= INTSORTSIZE) >= tmp_ptr; hi = lo)
+                  *hi = *lo;
+                *hi = c;
+              }
+          }
+      }
+  }
+}

integer_sort_mtypes.c

+/* This file is part of the Scotch distribution. It does
+** not have the stardard Scotch header with the INRIA
+** copyright notice because it is a very slight adaptation
+** of the qsort routine of glibc 2.4, taylored to match
+** Scotch needs. As Scotch is distributed according to the
+** CeCILL-C license, which is LGPL-compatible, no further
+** notices are required.
+*/
+
+/* Copyright (C) 1991,1992,1996,1997,1999,2004 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Written by Douglas C. Schmidt (schmidt@ics.uci.edu).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+/* If you consider tuning this algorithm, you should consult first:
+   Engineering a sort function; Jon Bentley and M. Douglas McIlroy;
+   Software - Practice and Experience; Vol. 23 (11), 1249-1265, 1993.  */
+
+
+/**
+ * @file integer_sort_mtypes.c
+ *
+ * File to include to create a multi-types sort function using qsort based
+ * algorithm. DO NOT compile directly.
+ *
+ * @author François Pellegrini
+ *
+ */
+#ifndef MAX_THRESH_2
+
+#define MAX_THRESH_2 6
+
+#define max_thresh_2                  (MAX_THRESH_2 * INTSORTSIZE(0)) /* Variable turned into constant */
+
+/* Stack node declarations used to store unfulfilled partition obligations. */
+typedef struct
+  {
+    char *lo;
+    char *hi;
+  } stack_node_2;
+
+/* The next 4 #defines implement a very fast in-line stack abstraction. */
+/* The stack needs log (total_elements) entries (we could even subtract
+   log(MAX_THRESH_2)).  Since total_elements has type size_t, we get as
+   upper bound for log (total_elements):
+   bits per unsigned char (CHAR_BIT) * sizeof(size_t).  */
+#define STACK_SIZE_2	(CHAR_BIT * sizeof (pastix_int_t))
+#define PUSH_2(low, high)	((void) ((top->lo = (low)), (top->hi = (high)), ++top))
+#define	POP_2(low, high)	((void) (--top, (low = top->lo), (high = top->hi)))
+#define	STACK_NOT_EMPTY_2	(stack < top)
+
+#endif /* MAX_THRESH_2 */
+
+/* Order size using quicksort.  This implementation incorporates
+   four optimizations discussed in Sedgewick:
+
+   1. Non-recursive, using an explicit stack of pointer that store the
+      next array partition to sort.  To save time, this maximum amount
+      of space required to store an array of SIZE_MAX is allocated on the
+      stack.  Assuming a 32-bit (64 bit) integer for size_t, this needs
+      only 32 * sizeof(stack_node) == 256 bytes (for 64 bit: 1024 bytes).
+      Pretty cheap, actually.
+
+   2. Chose the pivot element using a median-of-three decision tree.
+      This reduces the probability of selecting a bad pivot value and
+      eliminates certain extraneous comparisons.
+
+   3. Only quicksorts TOTAL_ELEMS / MAX_THRESH_2 partitions, leaving
+      insertion sort to order the MAX_THRESH_2 items within each partition.
+      This is a big win, since insertion sort is faster for small, mostly
+      sorted array segments.
+
+   4. The larger of the two sub-partitions is always pushed onto the
+      stack first, with the algorithm then concentrating on the
+      smaller partition.  This *guarantees* no more than log (total_elems)
+      stack size is needed (actually O(1) in this case)!  */
+
+/* To be defined :
+** INTSORTNAME : Name of function
+** INTSORTSIZE : Size of elements to sort
+** INTSORTSWAP : Swapping macro
+** INTSORTCMP  : Comparison function
+*/
+
+void
+INTSORTNAME (
+void ** const               pbase,                /*+ Array of arrays to sort   +*/
+const pastix_int_t                   total_elems)          /*+ Number of entries to sort +*/
+{
+  register char *base_ptr = (char *) (*pbase);
+
+  if (total_elems == 0)
+    /* Avoid lossage with unsigned arithmetic below.  */
+    return;
+
+  if (total_elems > MAX_THRESH_2)
+    {
+      char *lo = base_ptr;
+      char *hi = &lo[INTSORTSIZE(0) * (total_elems - 1)];
+      stack_node_2 stack[STACK_SIZE_2];
+      stack_node_2 *top = stack;
+
+      PUSH_2 (NULL, NULL);
+
+      while (STACK_NOT_EMPTY_2)
+        {
+          char *left_ptr;
+          char *right_ptr;
+
+          /* Select median value from among LO, MID, and HI. Rearrange
+             LO and HI so the three values are sorted. This lowers the
+             probability of picking a pathological pivot value and
+             skips a comparison for both the LEFT_PTR and RIGHT_PTR in
+             the while loops. */
+
+          char *mid = lo + INTSORTSIZE(0) * ((hi - lo) / INTSORTSIZE(0) >> 1);
+
+          if (INTSORTCMP (mid, lo))
+            INTSORTSWAP (mid, lo);
+          if (INTSORTCMP (hi, mid))
+            INTSORTSWAP (mid, hi);
+          else
+            goto jump_over;
+          if (INTSORTCMP ((void *) mid, (void *) lo))
+            INTSORTSWAP (mid, lo);
+        jump_over:;
+
+          left_ptr  = lo + INTSORTSIZE(0);
+          right_ptr = hi - INTSORTSIZE(0);
+
+          /* Here's the famous ``collapse the walls'' section of quicksort.
+             Gotta like those tight inner loops!  They are the main reason
+             that this algorithm runs much faster than others. */
+          do
+            {
+              while (INTSORTCMP ((void *) left_ptr, (void *) mid))
+                left_ptr += INTSORTSIZE(0);
+
+              while (INTSORTCMP ((void *) mid, (void *) right_ptr))
+                right_ptr -= INTSORTSIZE(0);
+
+              if (left_ptr < right_ptr)
+                {
+                  INTSORTSWAP (left_ptr, right_ptr);
+                  if (mid == left_ptr)
+                    mid = right_ptr;
+                  else if (mid == right_ptr)
+                    mid = left_ptr;
+                  left_ptr += INTSORTSIZE(0);
+                  right_ptr -= INTSORTSIZE(0);
+                }
+              else if (left_ptr == right_ptr)
+                {
+                  left_ptr += INTSORTSIZE(0);
+                  right_ptr -= INTSORTSIZE(0);
+                  break;
+                }
+            }
+          while (left_ptr <= right_ptr);
+
+          /* Set up pointers for next iteration.  First determine whether
+             left and right partitions are below the threshold size.  If so,
+             ignore one or both.  Otherwise, push the larger partition's
+             bounds on the stack and continue sorting the smaller one. */
+
+          if ((size_t) (right_ptr - lo) <= max_thresh_2)
+            {
+              if ((size_t) (hi - left_ptr) <= max_thresh_2)
+                /* Ignore both small partitions. */
+                POP_2 (lo, hi);
+              else
+                /* Ignore small left partition. */
+                lo = left_ptr;
+            }
+          else if ((size_t) (hi - left_ptr) <= max_thresh_2)
+            /* Ignore small right partition. */
+            hi = right_ptr;
+          else if ((right_ptr - lo) > (hi - left_ptr))
+            {
+              /* Push larger left partition indices. */
+              PUSH_2 (lo, right_ptr);
+              lo = left_ptr;
+            }
+          else
+            {
+              /* Push larger right partition indices. */
+              PUSH_2 (left_ptr, hi);
+              hi = right_ptr;
+            }
+        }
+    }
+
+  /* Once the BASE_PTR array is partially sorted by quicksort the rest
+     is completely sorted using insertion sort, since this is efficient
+     for partitions below MAX_THRESH_2 size. BASE_PTR points to the beginning
+     of the array to sort, and END_PTR points at the very last element in
+     the array (*not* one beyond it!). */
+
+#define min(x, y) ((x) < (y) ? (x) : (y))
+
+  {
+    char *const end_ptr = &base_ptr[INTSORTSIZE(0) * (total_elems - 1)];
+    char *tmp_ptr = base_ptr;
+    char *thresh = min(end_ptr, base_ptr + max_thresh_2);
+    register char *run_ptr;
+
+    /* Find smallest element in first threshold and place it at the
+       array's beginning.  This is the smallest array element,
+       and the operation speeds up insertion sort's inner loop. */
+
+    for (run_ptr = tmp_ptr + INTSORTSIZE(0); run_ptr <= thresh; run_ptr += INTSORTSIZE(0))
+      if (INTSORTCMP ((void *) run_ptr, (void *) tmp_ptr))
+        tmp_ptr = run_ptr;
+
+    if (tmp_ptr != base_ptr)
+      INTSORTSWAP (tmp_ptr, base_ptr);
+
+    /* Insertion sort, running from left-hand-side up to right-hand-side.  */
+
+    run_ptr = base_ptr + INTSORTSIZE(0);
+    while ((run_ptr += INTSORTSIZE(0)) <= end_ptr)
+      {
+        tmp_ptr = run_ptr - INTSORTSIZE(0);
+        while (INTSORTCMP ((void *) run_ptr, (void *) tmp_ptr))
+          tmp_ptr -= INTSORTSIZE(0);
+
+        tmp_ptr += INTSORTSIZE(0);
+        if (tmp_ptr != run_ptr)
+          {
+            char *trav;
+            int itertab;
+            trav = run_ptr + INTSORTSIZE(0);
+            while (--trav >= run_ptr)
+              {
+                char c = *trav;
+                char *hi, *lo;
+
+                for (hi = lo = trav; (lo -= INTSORTSIZE(0)) >= tmp_ptr; hi = lo)
+                  *hi = *lo;
+                *hi = c;
+              }
+
+            for (itertab = 1; itertab < INTSORTNTAB; itertab++)
+              {
+
+                trav = ((run_ptr - base_ptr)/INTSORTSIZE(0))
+                  *INTSORTSIZE(itertab) +
+                  (char *) pbase[itertab]+ INTSORTSIZE(itertab) ;
+
+                while (--trav >=  ((run_ptr - base_ptr)/INTSORTSIZE(0))
+                       *INTSORTSIZE(itertab) +
+                       (char *) pbase[itertab])
+                  {
+                    char c = *trav;
+                    char *hi, *lo;
+
+                    for (hi = lo = trav;
+                         (lo -= INTSORTSIZE(itertab)) >=
+                           ((tmp_ptr - base_ptr)/INTSORTSIZE(0))
+                           *INTSORTSIZE(itertab) +
+                           (char *) pbase[itertab]; hi = lo)
+                      *hi = *lo;
+                    *hi = c;
+                  }
+              }
+          }
+      }
+  }
+}