simple_dist_ijv_general_4x4_nospmcheck.c File Reference

Example demonstrating how to define a distributed sparse matrix out of a matrix allocated by the user in IJV format without requiring to call the expensive spmCheckAndCorrect function, and then solve a linear system with it. More...

Go to the source code of this file.

Detailed Description

Example demonstrating how to define a distributed sparse matrix out of a matrix allocated by the user in IJV format without requiring to call the expensive spmCheckAndCorrect function, and then solve a linear system with it.

Copyright

2025-2025 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria, Univ. Bordeaux. All rights reserved.

Version

6.4.0

Author

Florent Pruvost

Mathieu Faverge

Date

2025-09-02

   /
 
#include <mpi.h>
#include <pastix.h>
#include <spm.h>
#include <string.h>
/*
   Solve A*x = b, where
   A = [ 2 3 0 0 ]
       [ 0 4 5 0 ]
       [ 1 0 6 7 ]
       [ 0 0 8 9 ]
  
   b = [5, 9, 14, 17]
  
   solution must be
   x = [1, 1, 1, 1]
  
   This example reproduces the simple_dist_ijv_general_4x4 example with differences
   in the matrix initialization to remove the cost of the spmCheckAndCorrect function.
   First look at this previous example before examining this one.
  
   To avoid the call to spmCheckAndCorrect, it's important to meet the requirements
   of the PaStiX solver input:
     a) A CSC matrix (or CSR if distributed by row, see below)
   This can be easily done by a call to spmConvert if the matrix fulfills the
   following requirements:
     b) The rowptr (resp. colptr) is sorted for each column (resp. for each row)
     c) No duplicated entries exist
     d) If the matrix is general, the pattern is symmetrized by adding 0s (A+A^t)
  
   Thus we will store A as:
   A = [ 2 3 X 0 ]
       [ X 4 5 0 ]
       [ 1 X 6 7 ]
       [ 0 0 8 9 ]
   where X represents additional 0s that must be stored in the matrix to provide
   a symmetric pattern.
  
   Note that if the matrix is distributed by rows, you can use the SpmCSR format
   to store your problem and ask PaStiX to solve the transposed problem. This is
   shown in simple_trans.c
  
   The example is designed to work with 1 or 2 MPI processes.
  
  /
int main(int argc, char **argv)
{
    spmatrix_t    *spm;
    spm_int_t      nrhs = 1;
    int            rc   = 0;
    int            nmpi, rank;
    pastix_data_t *pastix_data = NULL; /*< Pointer to the storage structure required by PaStiX */
    pastix_int_t   iparm[IPARM_SIZE];  /*< Integer in/out parameters for PaStiX                */
    double         dparm[DPARM_SIZE];  /*< Floating in/out parameters for PaStiX               */
 
    MPI_Init( &argc, &argv );
    MPI_Comm_size( MPI_COMM_WORLD, &nmpi );
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    if ( (nmpi < 1) || (nmpi > 2) ) {
        fprintf(stderr, "This example is designed to run with 1 or 2 MPI processes.\n");
        return -1;
    }
 
    /**
       Initialize PaStiX parameters to default values
      /
    pastixInitParam( iparm, dparm );
 
    /**
       Get options from command line
       In this example, this is used only and only to read iparm/dparm values.
      /
    pastixGetOptions( argc, argv, iparm, dparm,
                      NULL, NULL, NULL, NULL );
 
    /**
       Initialize the sparse matrix
      /
    spm = malloc( sizeof( spmatrix_t ) );
    spmInitDist( spm, MPI_COMM_WORLD );
 
    spm->mtxtype = SpmGeneral; /* The matrix used is not symmetric or spd */
    spm->flttype = SpmDouble;
    spm->fmttype = SpmIJV;
    spm->baseval = 0;          /* C numbering is used */
    spm->dof     = 1;
 
    if ( nmpi > 1 ) {
        /*
           Make sure PaStiX/SpM knows if the matrix is replicated or distributed
           among the MPI processes.
           Note that with a single node, it is set to 1 by default.
          /
        spm->replicated = 0;
    }
 
    /*
       If the matrix is distributed among multiple nodes, loc2glob must be set
       to specify the distributions of the columns (or rows) among the nodes.
      /
    pastix_int_t *loc2glob = NULL;
    spm_int_t    *rowptr   = NULL;
    spm_int_t    *colptr   = NULL;
    double       *values   = NULL;
    double       *b        = NULL;
    double       *x        = NULL;
    if ( nmpi == 1 ) {
        /* If you run it with 1 MPI process, the matrix is fully defined in that process.
           A = [ 2 3 X 0 ]
               [ X 4 5 0 ]
               [ 1 X 6 7 ]
               [ 0 0 8 9 ]
          /
        spm->n   = 4;
        spm->nnz = 12;
        /*
           Note that these pointers (rowptr, colptr, values, dofs and loc2glob)
           could be allocated through a call to the spmAlloc function to make
           sure the allocations are compatible with free/alloc that may happen
           in the spmCheckAndCorrect. This is mandatory for C++ or Fortran code
           for example.
          /
        rowptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
        colptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
        values = (double *)   malloc(sizeof(double)    * (spm->nnz));
 
        /* Initialize A */
        rowptr[ 0] = 0; colptr[ 0] = 0; values[ 0] = 2.0; /* A[0,0] */
        rowptr[ 1] = 0; colptr[ 1] = 1; values[ 1] = 3.0; /* A[0,1] */
        rowptr[ 2] = 0; colptr[ 2] = 2; values[ 2] = 0.0; /* A[0,2] */
        rowptr[ 3] = 1; colptr[ 3] = 0; values[ 3] = 0.0; /* A[1,0] */
        rowptr[ 4] = 1; colptr[ 4] = 1; values[ 4] = 4.0; /* A[1,1] */
        rowptr[ 5] = 1; colptr[ 5] = 2; values[ 5] = 5.0; /* A[1,2] */
        rowptr[ 6] = 2; colptr[ 6] = 0; values[ 6] = 1.0; /* A[2,0] */
        rowptr[ 7] = 2; colptr[ 7] = 1; values[ 7] = 0.0; /* A[2,1] */
        rowptr[ 8] = 2; colptr[ 8] = 2; values[ 8] = 6.0; /* A[2,2] */
        rowptr[ 9] = 2; colptr[ 9] = 3; values[ 9] = 7.0; /* A[2,3] */
        rowptr[10] = 3; colptr[10] = 2; values[10] = 8.0; /* A[3,2] */
        rowptr[11] = 3; colptr[11] = 3; values[11] = 9.0; /* A[3,3] */
    }
    else if ( nmpi == 2 ) {
        /*
           If you run it with 2 MPI processes, the matrix is distributed among
           the two processes. PaStiX supports by default CSC matrices, thus IJV
           data distribution must be by columns (if distributed by row, you
           need to solve the transposed problem).
           Here, the first process has the first two columns, the second process
           has the last two columns. Note that any other column distribution
           could be used.
          /
        if ( rank == 0 ) {
            /* Local matrix on rank 0:
               col   0 1
               A = [ 2 3 ]
                   [ X 4 ]
                   [ 1 X ]
                   [ 0 0 ]
              /
            spm->n   = 2;
            spm->nnz = 6;
 
            loc2glob    = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->n));
            loc2glob[0] = 0; /* local index 0 corresponds to global index 0 */
            loc2glob[1] = 1; /* local index 1 corresponds to global index 1 */
 
            rowptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
            colptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
            values = (double *)   malloc(sizeof(double)    * (spm->nnz));
 
            /* Initialize A */
            rowptr[0] = 0; colptr[0] = 0; values[0] = 2.0; /* A[0,0] */
            rowptr[1] = 1; colptr[1] = 0; values[1] = 0.0; /* A[1,0] */
            rowptr[2] = 2; colptr[2] = 0; values[2] = 1.0; /* A[2,0] */
            rowptr[3] = 0; colptr[3] = 1; values[3] = 3.0; /* A[0,1] */
            rowptr[4] = 1; colptr[4] = 1; values[4] = 4.0; /* A[1,1] */
            rowptr[5] = 2; colptr[5] = 1; values[5] = 0.0; /* A[2,1] */
        }
        else {
            /* Local matrix on rank 1:
               col   2 3
               A = [ X 0 ]
                   [ 5 0 ]
                   [ 6 7 ]
                   [ 8 9 ]
              /
            spm->n   = 2;
            spm->nnz = 6;
 
            loc2glob    = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->n));
            loc2glob[0] = 2; /* local index 0 corresponds to global index 2 */
            loc2glob[1] = 3; /* local index 1 corresponds to global index 3 */
 
            rowptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
            colptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
            values = (double *)   malloc(sizeof(double)    * (spm->nnz));
 
            /* Initialize A */
            rowptr[0] = 0; colptr[0] = 2; values[0] = 0.0; /* A[0,2] */
            rowptr[1] = 1; colptr[1] = 2; values[1] = 5.0; /* A[1,2] */
            rowptr[2] = 2; colptr[2] = 2; values[2] = 6.0; /* A[2,2] */
            rowptr[3] = 3; colptr[3] = 2; values[3] = 8.0; /* A[3,2] */
            rowptr[4] = 2; colptr[4] = 3; values[4] = 7.0; /* A[2,3] */
            rowptr[5] = 3; colptr[5] = 3; values[5] = 9.0; /* A[3,3] */
        }
    }
    spm->rowptr   = rowptr;
    spm->colptr   = colptr;
    spm->values   = values;
    spm->loc2glob = loc2glob;
 
    /**
       Compute global N and NNZ
      /
    spmUpdateComputedFields( spm );
 
    /**
       Convert the spm to CSC format
      /
    rc = spmConvert( SpmCSC, spm );
    if ( rc != SPM_SUCCESS ) {
        spmExit( spm );
        fprintf( stderr, "Failed to convert the spm to CSC format\n" );
        return EXIT_FAILURE;
    }
 
    /* Print some properties */
    spmPrintInfo( spm, stdout );
    /* Print the matrix */
    spmPrint( spm, NULL );
 
    /*
       Now that the computed fields have been computed, we can create b and x
      
       x and b are of size nexp in case of multiple degree of freedom per
       unknown. In this example, nexp = n since dof = 1, and we could use
       any of the two fields.
      /
    x = (double *)malloc( sizeof(double) * spm->nexp );
    b = (double *)malloc( sizeof(double) * spm->nexp );
 
    if ( nmpi == 1 ) {
        /* Initialize b */
        b[0] =  5.0;
        b[1] =  9.0;
        b[2] = 14.0;
        b[3] = 17.0;
    }
    else {
        if ( rank == 0 ) {
            /* Initialize b */
            b[0] = 5.0;
            b[1] = 9.0;
        } else {
            /* Initialize b */
            b[0] = 14.0;
            b[1] = 17.0;
        }
    }
 
    /**
       Startup PaStiX
      /
    pastixInit( &pastix_data, MPI_COMM_WORLD, iparm, dparm );
 
    /**
       Perform ordering, symbolic factorization, and analyze steps
      /
    pastix_task_analyze( pastix_data, spm );
 
    /**
       Perform the numerical factorization
      /
    pastix_task_numfact( pastix_data, spm );
 
    /**
       Solve the linear system (and perform the optional refinement).
       PaStiX solve works with a single array for the right-hand side and solution.
       Let's copy b values in x before the solve.
      /
    memcpy( x, b, sizeof(double) * (spm->n) );
    pastix_task_solve( pastix_data, spm->nexp, nrhs, x, spm->nexp );
 
    /**
       Refine the solution if necessary.
      /
    if (iparm[IPARM_STATIC_PIVOTING] > 0 ) {
        fprintf( stdout,
                 " The number of static pivoting is %ld, which means that either the matrix is\n"
                 " either singular or close to singular, and static pivoting had to be used to\n"
                 " factorize the matrix. Iterative refinement must be applied to the solution\n"
                 " to refine it\n",
                 (long)iparm[IPARM_STATIC_PIVOTING] );
        /*
           Static pivoting is used, we can refine the solution to get better
           precision. The maximum number of iterations is controlled with
           iparm[IPARM_ITERMAX] (default = 250)
          /
        pastix_task_refine( pastix_data, spm->nexp, nrhs, b, spm->nexp, x, spm->nexp );
    }
 
    /**
       We can check the solution with the SpM interface. Warning: b stores Ax-b on exit
      /
    {
        double *bcpy = malloc( sizeof(double) * spm->nexp );
        memcpy( bcpy, b, sizeof(double) * spm->nexp );
        rc = spmCheckAxb( dparm[DPARM_EPSILON_REFINEMENT], nrhs, spm, NULL,
                          spm->nexp, bcpy, spm->nexp, x, spm->nexp );
        if ( rc != 0 ) {
            fprintf(stderr, "spmCheckAxb failed with error code: %d\n", rc);
            return rc;
        }
        free( bcpy );
    }
 
    /**
       Or we can directly compute b -= A * x, to check the correctness of the
       matrix, and check that the norm of the final b is null (at epsilon).
      /
    {
        rc = spmMatVec( SpmNoTrans, -1., spm, x, 1., b );
        if (rc != 0) {
            fprintf(stderr, "Matrix Vector multiplication failed with error code: %d\n", rc);
            return rc;
        }
        double resnorm = spmNormVec( SpmInfNorm, spm, b, 1 );
        fprintf( stdout, "The norm of || b - A * x || is %e\n", resnorm );
    }
 
    /* Destroy objects/contexts */
    pastixFinalize( &pastix_data );
 
    /* spmExit will free spm and also pointers rowptr, colptr, values, loc2glob */
    spmExit( spm );
    free( spm );
    free( x );
    free( b );
 
    MPI_Finalize();
 
    (void)argc;
    (void)argv;
    return rc;
}
/**

Definition in file simple_dist_ijv_general_4x4_nospmcheck.c.