PaStiX Handbook 6.4.0
Loading...
Searching...
No Matches
simple_dist_csc_general_4x4_reuse.c File Reference
PaStiX » Examples

Example demonstrating how to perform analysis on a common pattern matrix of multiple problems, and then solve multiple numerical problems re-using the same analysis. More...

Go to the source code of this file.

Detailed Description

Example demonstrating how to perform analysis on a common pattern matrix of multiple problems, and then solve multiple numerical problems re-using the same analysis.

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 A1 * x1 = b1, and A2 * x2 = b2, where
A1 = [ 2 3 0 0 ]
[ 0 4 5 0 ]
[ 1 0 6 7 ]
[ 0 0 8 9 ]
b1 = [5, 9, 14, 17]
A2 = [ 2 0 4 0 ]
[ 1 4 0 0 ]
[ 1 2 6 7 ]
[ 0 0 8 9 ]
b2 = [ 6, 5, 16, 17 ]
In both cases, the solution must be:
x1 = x2 = [1, 1, 1, 1]
Since the pattern of (A2 + A2^t) is included in the pattern of (A1 + A1^t),
the goal is to reuse as much computation as possible.
This example is inspired from the simple_dist_ijv_general_4x4_nospmcheck
example with the difference that we directly create the spm in the CSC format
to be able to reuse the spm.
WARNING 1: This is not possible if you use the IJV format that requires a
format change and thus a reallocation of some of the arrays.
WARNING 2: This is not possible if you use the spmCheckAndCorrect function
that may create the same issue as previously.
If you still want to do so, you'll have to traverse the prepared CSC to
update the values. So you need to be able to compute the values in your
matrix based on the coordinates (which is the case in many applications).
Thus, we will store A as the covering pattern of (A1 + A1^t) + (A2 + A2^t).
A = [ X X X 0 ]
[ X X X 0 ]
[ X X X X ]
[ 0 0 X X ]
where X represents additional zeros 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 = SpmCSC;
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 = [ X X X 0 ]
[ X X X 0 ]
[ X X X X ]
[ 0 0 X X ]
/
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.
/
colptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->n+1));
rowptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
values = (double *)malloc(sizeof(double) * (spm->nnz));
/* Initialize A */
colptr[0] = 0;
colptr[1] = 3;
colptr[2] = 6;
colptr[3] = 10;
colptr[4] = 12;
/* Initialize A1 */
rowptr[ 0] = 0; values[ 0] = 2.0; /* A[0,0] */
rowptr[ 1] = 1; values[ 1] = 0.0; /* A[1,0] */
rowptr[ 2] = 2; values[ 2] = 1.0; /* A[2,0] */
rowptr[ 3] = 0; values[ 3] = 3.0; /* A[0,1] */
rowptr[ 4] = 1; values[ 4] = 4.0; /* A[1,1] */
rowptr[ 5] = 2; values[ 5] = 0.0; /* A[2,1] */
rowptr[ 6] = 0; values[ 6] = 0.0; /* A[0,2] */
rowptr[ 7] = 1; values[ 7] = 5.0; /* A[1,2] */
rowptr[ 8] = 2; values[ 8] = 6.0; /* A[2,2] */
rowptr[ 9] = 3; values[ 9] = 8.0; /* A[3,2] */
rowptr[10] = 2; values[10] = 7.0; /* A[2,3] */
rowptr[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 only supports distributed matrices in CSC
format such that the distribution must be by columns, meaning that all
elements in a column must be present in the same process. Here, the
first process has the first two columns, the second process has the
last two columns. Cyclic distribution could also be used here,
alternating columns between the two processes.
/
if ( rank == 0 ) {
/* Local matrix on rank 0:
col 0 1
A = [ X X ]
[ X X ]
[ X 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 */
colptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->n+1));
rowptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
values = (double *)malloc(sizeof(double) * (spm->nnz));
/* Initialize A */
colptr[0] = 0;
colptr[1] = 3;
colptr[2] = 6;
/* Initialize A1 */
rowptr[0] = 0; values[0] = 2.0; /* A[0,0] */
rowptr[1] = 1; values[1] = 0.0; /* A[1,0] */
rowptr[2] = 2; values[2] = 1.0; /* A[2,0] */
rowptr[3] = 0; values[3] = 3.0; /* A[0,1] */
rowptr[4] = 1; values[4] = 4.0; /* A[1,1] */
rowptr[5] = 2; values[5] = 0.0; /* A[2,1] */
}
else {
/* Local matrix on rank 1:
col 2 3
A = [ X 0 ]
[ X 0 ]
[ X X ]
[ X X ]
/
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 */
colptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->n+1));
rowptr = (spm_int_t *)malloc(sizeof(spm_int_t) * (spm->nnz));
values = (double *)malloc(sizeof(double) * (spm->nnz));
/* Initialize A */
colptr[0] = 0;
colptr[1] = 4;
colptr[2] = 6;
/* Initialize A1 */
rowptr[0] = 0; values[0] = 0.0; /* A[0,2] */
rowptr[1] = 1; values[1] = 5.0; /* A[1,2] */
rowptr[2] = 2; values[2] = 6.0; /* A[2,2] */
rowptr[3] = 3; values[3] = 8.0; /* A[3,2] */
rowptr[4] = 2; values[4] = 7.0; /* A[2,3] */
rowptr[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 );
/* Print some properties */
spmPrintInfo( spm, stdout );
if (rank == 0) {
fprintf(stdout,
"+-------------------------------------------------+\n"
"+ Solve the first problem +\n"
"+-------------------------------------------------+\n");
}
/* 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 b1 */
b[0] = 5.0;
b[1] = 9.0;
b[2] = 14.0;
b[3] = 17.0;
}
else {
if ( rank == 0 ) {
/* Initialize b1 */
b[0] = 5.0;
b[1] = 9.0;
} else {
/* Initialize b1 */
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 );
}
if ( nmpi == 1 ) {
fprintf(stdout, "The solution found is [%e, %e, %e, %e]\n",
x[0], x[1], x[2], x[3]);
}
else {
if ( rank == 0 ) {
fprintf(stdout, "The solution found is [%e, %e, ",
x[0], x[1]);
}
MPI_Barrier(MPI_COMM_WORLD);
if ( rank == 1 ) {
fprintf(stdout, "%e, %e]\n",
x[0], x[1]);
}
}
/**
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 );
}
if (rank == 0) {
fprintf(stdout,
"+-------------------------------------------------+\n"
"+ Solve the second problem +\n"
"+-------------------------------------------------+\n");
}
/*
Solve the second problem
/
if (nmpi == 1) {
/* Initialize A2 */
/* rowptr[ 0] = 0; */ values[ 0] = 2.0; /* A[0,0] */
/* rowptr[ 1] = 1; */ values[ 1] = 1.0; /* A[1,0] */
/* rowptr[ 2] = 2; */ values[ 2] = 1.0; /* A[2,0] */
/* rowptr[ 3] = 0; */ values[ 3] = 0.0; /* A[0,1] */
/* rowptr[ 4] = 1; */ values[ 4] = 4.0; /* A[1,1] */
/* rowptr[ 5] = 2; */ values[ 5] = 2.0; /* A[2,1] */
/* rowptr[ 6] = 0; */ values[ 6] = 4.0; /* A[0,2] */
/* rowptr[ 7] = 1; */ values[ 7] = 0.0; /* A[1,2] */
/* rowptr[ 8] = 2; */ values[ 8] = 6.0; /* A[2,2] */
/* rowptr[ 9] = 3; */ values[ 9] = 8.0; /* A[3,2] */
/* rowptr[10] = 2; */ values[10] = 7.0; /* A[2,3] */
/* rowptr[11] = 3; */ values[11] = 9.0; /* A[3,3] */
/* Initialize b2 */
b[0] = 6.0;
b[1] = 5.0;
b[2] = 16.0;
b[3] = 17.0;
}
else if (nmpi == 2) {
if (rank == 0) {
/* Initialize A2 */
/* rowptr[0] = 0; */ values[0] = 2.0; /* A[0,0] */
/* rowptr[1] = 1; */ values[1] = 1.0; /* A[1,0] */
/* rowptr[2] = 2; */ values[2] = 1.0; /* A[2,0] */
/* rowptr[3] = 0; */ values[3] = 0.0; /* A[0,1] */
/* rowptr[4] = 1; */ values[4] = 4.0; /* A[1,1] */
/* rowptr[5] = 2; */ values[5] = 2.0; /* A[2,1] */
/* Initialize b2 */
b[0] = 6.0;
b[1] = 5.0;
}
else {
/* Initialize A2 */
/* rowptr[0] = 0; */ values[0] = 4.0; /* A[0,2] */
/* rowptr[1] = 1; */ values[1] = 0.0; /* A[1,2] */
/* rowptr[2] = 2; */ values[2] = 6.0; /* A[2,2] */
/* rowptr[3] = 3; */ values[3] = 8.0; /* A[3,2] */
/* rowptr[4] = 2; */ values[4] = 7.0; /* A[2,3] */
/* rowptr[5] = 3; */ values[5] = 9.0; /* A[3,3] */
/* Initialize b2 */
b[0] = 16.0;
b[1] = 17.0;
}
}
/* Print the matrix */
spmPrint(spm, NULL);
/**
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 );
}
if ( nmpi == 1 ) {
fprintf(stdout, "The solution found is [%e, %e, %e, %e]\n",
x[0], x[1], x[2], x[3]);
}
else {
if ( rank == 0 ) {
fprintf(stdout, "The solution found is [%e, %e, ",
x[0], x[1]);
}
MPI_Barrier(MPI_COMM_WORLD);
if ( rank == 1 ) {
fprintf(stdout, "%e, %e]\n",
x[0], x[1]);
}
}
/**
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;
}
/**
pastix_int_t
BEGIN_C_DECLS typedef int pastix_int_t
Definition datatypes.h:51
pastixFinalize
void pastixFinalize(pastix_data_t **pastix_data)
Finalize the solver instance.
Definition api.c:928

Definition in file simple_dist_csc_general_4x4_reuse.c.