Mat: Matrices in PETSc¶
This tutorial is an introduction to matrix objects (Mat) in PETSc.
Working with a matrix in PETSc involves several steps:
Create a matrix
Insert or add values into the matrix
Process (or assemble) the matrix
Use the matrix for computations
Destroy the matrix
In this tutorial, we demonstrate how to create a tridiagonal parallel matrix. A matrix of this form can arise, for example, from a three point finite difference discretization of the one-dimensional Laplace equation. We will also show basic matrix and matrix-vector operations. An example of creating a sequential tridiagonal matrix can be found in PETSc KSP tutorial ex1 and a parallel example is shown in PETSc KSP tutorial ex23.
Creating a PETSc matrix¶
We start by setting up a PETSc program.
To utilize the matrix class, we need to include petscmat.h which
automatically includes header files for basic PETSc routines, index sets,
and viewers.
We also need to declare variables to hold the data we will use in this tutorial.
In this case, we declare a Mat and two Vec objects, along with other
PetscInt and PetscScalar variables that will be used to, for example, set matrix values.
PETSc programs usually begin with PetscInitialize, which sets up MPI, among other things.
static char help[] = "Create a tridiagonal matrix.\n\n";
#include <petscmat.h>
int main(int argc,char **args)
{
Mat A;
Vec x,b;
PetscInt i,n=10,col[3],rstart,rend,nlocal,mlocal;
PetscErrorCode ierr;
PetscScalar value[3];
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
We create a matrix with the MatCreate function.
This will create a sequential matrix if the program is executed with one processor
or a parallel matrix if it is executed with more than one processor.
PETSc supports a variety of matrix implementations since no single matrix format is appropriate for all problems.
For example, sparse matrices can be stored in compressed row (AIJ) format.
The matrix type can be coded into a program or set from the command line using, for example, -mat_type mpiaij
in the same way that vector formats can be set.
We create an n by n square matrix, let PETSc decide how to partition it across processors,
and check for command line options.
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
To insert or add entries to a matrix, we use the MatSetValues function.
A processor can set any entry in a matrix even if the entry is not locally owned,
but for efficient assembly it is important that each processor
contributes to entries it owns locally (as much as is possible).
To accomplish this, we identify the portion of the matrix on each processor
and insert values into the matrix.
In this example, we insert values one row at a time in a loop by providing the column indices and values.
We handle the boundaries separately since these rows only contain two entries.
The final step before we can use the matrix is assembly.
ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
/* Boundary */
if (!rstart) {
rstart = 1;
i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
if (rend == n) {
rend = n-1;
i = n-1; col[0] = n-2; col[1] = n-1; value[0] = -1.0; value[1] = 2.0;
ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
/* Interior */
value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
for (i=rstart; i<rend; i++) {
col[0] = i-1; col[1] = i; col[2] = i+1;
ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
Creating a compatible vector¶
Matrices are often used with vectors in computations. When we create a parallel vector to be used with a matrix, we can force it to have the same layout by specifying the local and global size of the vector. Note that we can create one vector from scratch and duplicate it as needed.
ierr = MatGetLocalSize(A,&nlocal,&mlocal);CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecSetSizes(x,mlocal,n);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
ierr = VecSet(x,1.0);CHKERRQ(ierr);
Matrix operations in PETSc¶
PETSc supports many matrix routines such as those that compute the trace of a matrix and matrix norms. There are also routines that operate on matrices like scaling a matrix or computing its conjugate. Many matrix-vector operations are supported, as well. A full list of matrix routines can be found here and a few are shown below.
ierr = MatMult(A,x,b);CHKERRQ(ierr);
ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = MatScale(A,2.0);CHKERRQ(ierr);
ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
All PETSc objects should be destroyed when they are no longer needed, so we conclude by freeing work space and finalizing the program.
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
Running this program with two MPI ranks produces the output below.
Vec Object: 2 MPI processes
type: mpi
Process [0]
1.
0.
0.
0.
0.
Process [1]
0.
0.
0.
0.
1.
Mat Object: 2 MPI processes
type: mpiaij
row 0: (0, 4.) (1, -2.)
row 1: (0, -2.) (1, 4.) (2, -2.)
row 2: (1, -2.) (2, 4.) (3, -2.)
row 3: (2, -2.) (3, 4.) (4, -2.)
row 4: (3, -2.) (4, 4.) (5, -2.)
row 5: (4, -2.) (5, 4.) (6, -2.)
row 6: (5, -2.) (6, 4.) (7, -2.)
row 7: (6, -2.) (7, 4.) (8, -2.)
row 8: (7, -2.) (8, 4.) (9, -2.)
row 9: (8, -2.) (9, 4.)