simple_dist_csr_symmetric_lower_4x4.c File Reference

Example demonstrating how to define a distributed sparse symmetric matrix from a matrix allocated by the user in CSR format and 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 symmetric matrix from a matrix allocated by the user in CSR format and 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

Date

2025-09-02

   /
 
#include <mpi.h>
#include <pastix.h>
#include <spm.h>
#include <string.h>
/*
   Solve A*x = b, where
   A = [ 2 . . . ]
       [ 1 3 . . ]
       [ 0 1 4 . ]
       [ 1 0 1 2 ]
  
   b = [4, 5, 6, 4]
  
   solution must be
   x = [1, 1, 1, 1]
  
   We demonstrate how to enter the symmetric matrix using the CSR format.
  
   The matrix must be defined in a way that the distribution is by rows
   (i.e., all elements in a row must be present in the same process).
   PaStiX only supports distributed matrices in CSC format internally, so we
   need to convert the CSR matrix to CSC format.
  
   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 = SpmSymmetric; /* The matrix used is symmetric */
    spm->flttype = SpmDouble;
    spm->fmttype = SpmCSR;
    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 . . . ]
               [ 1 3 . . ]
               [ 0 1 4 . ]
               [ 1 0 1 2 ]
          /
        spm->n   = 4;
        spm->nnz = 8;
        /*
           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->n + 1));
        colptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
        values = (double *)   malloc(sizeof(double)    * (spm->nnz));
 
        /*
           Initialize A:
           We give only the lower part as the matrix is symmetric.
           It works also if the upper part is given.
          /
        rowptr[0] = 0; rowptr[1] = 1; rowptr[2] = 3; rowptr[3] = 5; rowptr[4] = 8;
        colptr[0] = 0; values[0] = 2.0; /* A[0,0] */
        colptr[1] = 0; values[1] = 1.0; /* A[1,0] */
        colptr[2] = 1; values[2] = 3.0; /* A[1,1] */
        colptr[3] = 1; values[3] = 1.0; /* A[2,1] */
        colptr[4] = 2; values[4] = 4.0; /* A[2,2] */
        colptr[5] = 0; values[5] = 1.0; /* A[3,0] */
        colptr[6] = 2; values[6] = 1.0; /* A[3,2] */
        colptr[7] = 3; values[7] = 2.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 the
           data distribution must be by columns, meaning that all elements in a
           column must be present in the same process. Here, by symmetry of the
           CSR format, the first process has the first two rows, and the second
           process has the last two rows.
          /
        if ( rank == 0 ) {
            /* Local matrix on rank 0:
                               row
               A = [ 2 . . . ] 0
                   [ 1 3 . . ] 1
              /
            spm->n   = 2;
            spm->nnz = 3;
 
            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 */
 
            colptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
            rowptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->n + 1));
            values = (double *)   malloc(sizeof(double)    * (spm->nnz));
 
            /*
               Initialize A:
               We give only the lower part as the matrix is symmetric.
              /
            rowptr[0] = 0; rowptr[1] = 1; rowptr[2] = 3;
            colptr[0] = 0; values[0] = 2.0; /* A[0,0] */
            colptr[1] = 0; values[1] = 1.0; /* A[1,0] */
            colptr[2] = 1; values[2] = 3.0; /* A[1,1] */
        }
        else {
            /* Local matrix on rank 1:
                               row
               A = [ 0 1 4 . ] 2
                   [ 1 0 1 2 ] 3
              /
            spm->n   = 2;
            spm->nnz = 5;
 
            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->n + 1));
            colptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
            values = (double *)   malloc(sizeof(double)    * (spm->nnz));
 
            /*
               Initialize A:
               We give only the lower part as the matrix is symmetric.
              /
            rowptr[0] = 0; rowptr[1] = 2; rowptr[2] = 5;
            colptr[0] = 1; values[0] = 1.0; /* A[2,1] */
            colptr[1] = 2; values[1] = 4.0; /* A[2,2] */
            colptr[2] = 0; values[2] = 1.0; /* A[3,0] */
            colptr[3] = 2; values[3] = 1.0; /* A[3,2] */
            colptr[4] = 3; values[4] = 2.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] = 4.0;
        b[1] = 5.0;
        b[2] = 6.0;
        b[3] = 4.0;
    }
    else {
        if ( rank == 0 ) {
            /* Initialize b */
            b[0] = 4.0;
            b[1] = 5.0;
        } else {
            /* Initialize b */
            b[0] = 6.0;
            b[1] = 4.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_csr_symmetric_lower_4x4.c.